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1 . INTRODUCTION 


SAMCEP  is  a Monte  Carlo  neutron  and  secondary  gamma  transport  code  which 
is  designed  to  calculate  differences  in  fluxes  and  in  flux-derived  quantities 
such  as  dose,  as  functions  of  differences  in  nuclear  data.  Source  spectral 
changes  can  also  be  considered.  SAMCEP  is  an  extension  of  SAM-CE1 , a general 

purpose  three-dimensional  Monte  Carlo  transport  code.  SAMCEP,  originally  de- 

★ 

veloped  under  contract  DAAD-05-70-C-0295  , uses  the  latest  (ENDF/B)  cross  sec- 
2 4 

tions  ' to  a high  degree  of  precision 

The  extension  of  SAM-CE  to  SAMCEP  involves  the  addition  of  the  capability 
of  running  several  similar  Monte  Carlo  problems  simultaneously.  The  mathematical 
basis  of  this  extension  is  the  simultaneous  use,  in  several  correlated  problems, 
of  the  same  Monte  Carlo  histories,  taking  account  of  the  effects  of  the  cross 
section  differences  via  a series  of  problem-dependent  weights  associated  with 
each  history. 

Using  SAMCEP,  flux  differences  due  to  small  cross  section  variations  can 
be  calculated  much  more  rapidly  and  accurately  than  by  conventional  methods, 
i.e.,  by  independent  computations.  The  process  of  measuring,  reducing,  evalua- 
ting, and  communicating  nuclear  data  is  slow  and  expensive.  It  is  essential, 
therefore,  for  proper  management  of  such  an  effort  to  know  whether  errors  or  un- 
certainties in  data  result  in  operationally  significant  errors  or  uncertainties 
in  transport  computations.  SAMCEP  is  a practical  tool  that  can  provide  inform- 
ation permitting  one  to  direct  cross  section  research  to  obtaining  more  accurate 
data  in  those  areas  where  the  lack  of  precision  leads  to  the  greatest  uncertainty 
in  the  transport  calculation.  ■' 


* 

The  bounded  f lux-at-a-point  capability  was  developed  under  DAAD-05-73-C-0072 
The  secondary  gamma  capability  under  DAAD-05-75-C-0735. 
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In  order  to  describe  the  capability  of  SAMCEP,  it  is  necessary  first  to 
define  a perturbation:  a perturbation  is  a specific  change  in  some  portion  of 

the  neutron  source  or  changes  in  the  nuclear  data  for  any  or  all  of  the  ele- 
ments that  make  up  a medium  in  a transport  problem,,  SAMCEP  can  handle,  in  a 

single  calculation,  up  to  ten  correlated  problems,  each  correlated  problem  in- 

* 

eluding  a combination  of  up  to  ten  perturbations  . The  number  of  allowed  per- 
turbations of  the  neutron  data  is,  therefore,  one  hundred,  (Perturbations  of 


the  gamma  production  data  are  described  in  the  chapter  on  PROGRAM  SAMGAM) . 

Perturbations  of  the  neutron  data  are  classified  by  type  as  follows: 

Type  1 composition  (i.a„,  concentration  perturbatron)  ,- 
these  perturbations  carry  over  to  a secondary 
gamma  problem; 

Type  2 the  complete  set  of  cross  section  data  of  an 
element  (i„e.,  an  alternate  set  of  cross 
section  data  of  a specific  element) ; 

Type  3 change  of  microscopic  total  cross  sections; 

Type  4 change  of  microscopic  scattering  cross 

sections; 

Type  5 change  of  microscopic  inelastic  scattering 
cross  sections; 

Type  6 change  in  angular  distribution  of  elastic 
scattering ; 

Type  7 change  in  secondary  energy  distribution  in 
continuum  inelastic  scattering; 

Type  8 change  in  cross  sections  for  inelastic  level 
excitation; 

Type  9 change  in  angular  distribution  of  discrete 
level  inelastic  scattering; 

Type  10  change  in  source  spectrum 
Types  3-9  are  defined  within  specific  energy  ranges. 


exclusive  of  perturbations  of  the  gamma  production  data  for  a secondaiy 
gamma  problem 
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SAMCEP  consists  of  five  separate  programs  which  may  be  executed  as  separ- 
ate jobs  or  in  sequence  as  one  job.  Information  is  passed  from  one  program  to 
another  either  by  physical  tapes  or  by  disk  files.  (Detailed  descriptions  of 
each  program  and  its  input  data  requirements  will  be  given  below  in  separate 
chapters.)  The  five  programs  of  SAMCEP  are  (in  the  order  in  which  they  are 
executed) : 

1.  Program  SAMIN  (primary  neutron  problem)  - This  program  reads 

the  input  data  describing  the  neutron  perturbations  and 

converts  them  into  the  same  format  as  the  output  of  program 

■f 

SAM-X  (processor  of  ENDF/B  cross  section  data  ) . 

2.  Program  SAMSAM  (primary  neutron  problem)  - This  program  reads 

in  the  data  generated  by  SAMIN  from  a tape  or  a disk  file, 

* * * 

generates  sampling  CHI-  and  ENN-  tables  when  such  perturba- 
tions exist,  (type  2,  6 and  7) , and  writes  such  information 
on  a tape  or  disk  file.  The  other  basic  function  of  SAMSAM 
is  to  process  cross  section  data  into  energy  band  structure 
by  calling  subroutine  BAND. 


The  processor  SAM-X  is  a component  of  the  SAM-CE  system  of  programs 

(Ref.  1) 


Table  of  CHI-boundaries,  CHI= (l-cos0)/2,  where  $ is  scattering  angle  of 
equi-probable  CHI-bins  following  elastic  scattering. 

Table  of  energy  boundaries  for  equi-probabl  •jmerging-energy  bins  in 
continuum  inelastic  scattering. 
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3. 


Program  SAMCAR  (primary  neutron  and  secondary  gamma) 


This  is  the  main  Monte  Carlo  program.  By  calls  to  a series 
of  independent  routines,  it  performs  the  following  basic 
functions: 

a)  process  geometry  data  by  calling  subroutine  GENI, 

b)  process  the  unperturbed  and  perturbed  (if  any) 
source  spectrum  by  calling  subroutine  SOUCAL, 

c)  perform  the  transport  calculations  and  score  the 
fluxes  by  calling  subroutine  CARLO. 

★ 

SAMCAR  also  writes  the  fluxes  by  supergroups  of  all  problems  (unperturbed 
and  perturbed)  fc  r-.c'i  aggregate  on  a tape  or  a disk  file  to  be  edited  by  the 
next  program  SAMOUT. 

All  neutron  perturbation  data  will  be  stored  in  core  memory  all  the  time 
(except  for  perturbation  type  2)  and  in  the  same  format  es  the  unperturbed  data 
which  is  the  output  of  SAM-X  (processor  of  ENDF/B  data)..  For  perturbation  type  2, 
where  an  alternate  set  of  cross  section  data  replaces  the  unperturbed  set,  the 
input  is  required  to  have  the  ENDF/B  format.  SAM-X  and  BAND  will  process  this  as 
an  additional  element.  (SAMCAR  treats  the  data  as  an  alternate  set  for  the  same 
element. ) 

For  secondary  gamma  problems,  perturbations  of  the  gamma  production  data 
(if  any)  are  incorporated  into  the  external  source  tape  generated  by  program 
SAMGAM. 


The  term  supergroup  refers  to  a subdivision  of  the  output  energy  range, 
analogous  to  a band , which  refers  to  a cross  section  energy  subdivision 
Particles  are  tracked  within  a superbin,  defined  as  the  intersection  of 
the  current  band  and  current  supergroup.  (Refer  to  Sec.  5.2.6) 
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4.  Program  SAMOUT  (primary  neutron  and  secondary  gamma) 


I 


This  program  reads  in  the  aggregate  tape  and  processes  a 
statistical  tape.  It  then  edits  the  results  of  all  prob- 
lems, i.e.,  it  displays  these  results  and  deviation  estimates 
of  all  the  problems,  and  also  the  differences  of  these  re- 
sults and  the  deviation  estimates  of  their  differences  for 
all  specified  pairs  of  problems 

5 Program  SAMGAM  (secondary  gamma  problem)  - This  is  the  pre- 
processor for  a secondary  gamma  problem.  It  presumes  the 
prior  execution  of  the  (SAMIN,  SAMSAM,  SAMCAR,  SAMOUT)  se- 
quence for  the  precursor  primary  neutron  problem,  which  gener- 
ates a tape  of  non-elastic  (gamma-producing)  neutron  inter- 
actions Its  principal  functions,  performed  by  calling  several 
subroutines,  are: 

a)  alter  neutron  perturbation  tape,  retaining  only  con- 
centration perturbations  (if  any),  by  calling  PREPS; 

b)  pr  pare  organized  gamma  interaction  cross  section 
tape  by  calling  BANDG; 

c)  read  in  user  prepared  perturbations  of  the  gamma  pro- 
duction data,  and 

d)  generate  external  source  tape,  utilizing  the  interaction 
tape  of  the  precursor  neutron  primary  and  the  (possibly 
perturbed)  gamma  production  data,  by  calling  SAMSOU. 

The  execution  of  program  SAMGAM  is  followed  by  executions  of  SAMCAR  (in 
1 ho  secondary  gamma  mode)  and  SAMOUT  to  complete  the  secondary  gamma  sequence. 
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2.  TECHNICAL  BACKGROUND 


2.1  Simultaneous  Tracking 

In  order  to  estimate  accurately  the  differences  in  problem  results,  each 
Monte  Carlo  history  is  used  to  estimate  all  problems  simultaneously.  As  a 
result  there  is  a high  degree  of  correlation  between  the  answers  computed  for 
the  individual  problems.  Hence,  the  percentage  fluctuations  in  the  differences, 
due  to  the  randomness  inherent  in  the  Monte  Carlo,  are  comparable  to  those  in 
the  individual  answers.  Therefore,  the  absolute  error  in  the  estimated  differ- 
ence is  significantly  smaller  than  the  error  in  the  difference  of  the  results 
of  separate,  uncorrelated  calculations. 

The  implementation  of  this  procedure,  the  use  of  the  same  history  for 
several  problems,  sometimes  called  correlated  sampling,  is  based  on  the  fact 
that  it  is  always  possible  to  carry  out  a Monte  Carlo  calculation  using  prob- 
ability distributions  other  than  those  that  describe  the  natural  stochastic 
processes  of  radiation.  In  doing  so,  a cumulative  product  of  ratios  of  natural 
to  altered  distributions  is  calculated  as  a weight  with  which  an  answer  is 
counted.  Thus,  it  is  possible  to  use  as  'natural'  and  'altered'  distributions 
those  derived  from  alternative  data.  Several  answers  may  be  obtained;  one  for 
the  distribution  actually  sampled  (with  weight  equal  to  one)  and  others  from  the 
ratios  of  distributions.  This  can  be  generalized  by  sampling  from  probability 
density  functions  not  corresponding  to  any  physical  problem.  The  sampling 
density  can  be  defined  as  the  one  which  attempts  to  minimize  the  statistical  un- 
certainties in  the  difference  of  the  results  of  the  correlated  problems. 

In  SAMCEP,  the  sampling  transport  distribution  is  defined  as  the  least 
upper  bound  of  the  individual  transport  distributions.  All  other  sampling 
distributions  (source  sampling,  and  sampling  of  interaction  events)  are  proper 
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probability  distribution  functions,,  They  are,  therefore,  defined  as  the 

* 

renormalized  least  upper  bound  (l.u.b„)  of  the  individual  distributions  . 

The  various  sampling  distributions  are  described  in  more  detail  below. 

Transport  Sampling  Distribution 

In  any  region,  the  natural  transport  kernel  for  the  ith  problem  is  of 
the  form  Xexp (-p 1s)ds,  where  p 1 is  the  total  cross  section  at  the  current 

energy  E,  for  concentrations  and  cross  sections  appropriate  for  the  ith  problem. 
The  F^'s  reflect  the  previous  weight  adjustments  and  transport  through  pre- 
vious regions* 

The  ordinary  sampling  distribution  for  all  problems  i,  i=l,n,  is  the 
least  upper  bound  or  envelope  of  the  individual  distributions  Fp  1exp (-!■  * s) d , 
which  is  piece-wise  exponential  and  can,  therefore,  be  selected  by  the  tundard 
methods*  For  BFAP  estimation,  this  procedure  is  modified  (see  Appendix  E) . 

* * 

Source  Sampling  Distribution  (neutron  transport) 

Each  problem  i,  i=l,n,  can  have  a perturbed  source  energy  distribution 
S^(E).  The  sampling  distribution  is  constructed  by  first  construe:  irg  • 
which  is  the  envelope  of  all  S^(E),  biasing  it  by  constructing  S(E)  v.^  i wl  (o 
W^,  (E)  are  the  energy-dependent  weights  in  the  source  region,  and  fine 
normalizing,  obtaining  the  sampling  distribution 


S(E) 


S(E)/W  (E) 

fci 


(S  (E’  ) /W  (E ' ) ) dE ' 

Cj 


The  use  of  the  least  upper  bound  of  the  distributions,  called  the  envelope,  is 
motivated  by  the  following  considerations:  it  minimizes  the  maximum  relative 
weight  factor  in  sampling,  and  so  tends  to  prevent  large  variances;  and  it  can 
be  constructed  automatically  in  the  course  of  the  calculation,  so  that  a de- 
tailed examination  of  the  perturbed  data  is  not  necessary. 

**An  external  source  tape  is  generated  by  SAMGAM  for  secondary  gamma  transport . 
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Interaction  Sampling  Distributions 

Selection  of  Element:  Given  an  interaction  in  a region,  the  first  samp- 

ling implemented  is  that  of  the  particular  element  with  which  the  collision 
occurred 
Let 

li1  = macroscopic  total  cross  section, 
i 

o.  = microscopic  total  cross  section  for  element  j, 

C1  = element  concentration  of  jth  element, 

3 

NE  = total  number  of  elements 


NE  . . 

U = £ C .o . , for  all  i (i=l,  total  number  of  problems) 
j=l  3 3 


Let  p^  be  the  probability  of  picking  element  j for  the  i*"*1  problem: 

p1  = -JL-1 

3 U1  • 


Then  the  sampling  probability  p^  for  picking  element  j for  collision  is  the  re- 


normalized 1 u.b.  of  the  p. 

3 


F.C^o1 

pS  = max.  JLJLJ. 

3 i ! 


NORM  for  all  j. 


where  NORM  = £ MAX. 


l i 

! 

. 1,1  J 
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Selection  of  Interaction  Type;  Given  an  interaction  with  a particular  element. 


the  next  selection  is  the  type  of  interaction  to  occur.  If  the  microscopic 

elemental  total,  scattering  and  inelastic  cross  section  data  are  given  by 

ot , o1  and  cr*  , then 
t'  s in' 


where 


1 

o 

a 


i 

CT 


S 


are  the  probabilities  of  having  absorption,  elastic  scattering,  and  inelastic 
scattering,  respectively,  for  the  i*”*1  problem.  We  denote  the  individual  distri- 
bution (the  ith  distribution)  by 

i 

a 

p = — r , where  r may  denote  absorption,  elastic  or  inelastic 
r ot  scattering  reactions. 


Then,  again  by  choosing  the  envelope  of  this  discrete  function,  we  have  for  the 
sampling  probability  distribution  for  picking  reaction  r 

MAX. iF.p1] 
i i*r 

s 

Pr  = r T~ 

1 MAXi(FiPrl 

th  i s 

The  perturbation  weight  F^  for  the  i problem  is  then  multiplied  by  p /p‘ 

Since  we  allow  perturbation  of  the  c^n  cross  section  and  also  the  cross 
section  for  level  excitation  (type  2,  5 and  8),  we  have  to  compute  the  sampling 
distribution  for  picking  a level  or  continuum.  This  sampling  distribution  is 
again  the  envelope  of  all  the  individual  distributions.  When  a certain  level 
or  continuum  is  chosen,  the  weight  is  adjusted  in  the  same  manner  as  described 
above . 
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Given  an  elastic  scattering,  the  sampling  distribution  for  the  angle  of 
scattering  is  again  constructed  as  the  envelope  of  the  angular  distributions 
corresponding  to  the  different  perturbations,  and  is  renormalized.  The  same 
applies  to  the  energy  distribution  of  inelastic  scattering  leading  to  contin- 
uous spectra, 

2,2  On  Energy  Binning  of  Perturbations 

In  the  course  of  the  main  Monte  Carlo,  at  every  stage  at  which  a decision 
is  to  be  made  (e„g„,  distance  to  a collision  element  with  which  a collision 
occurs,  reaction  type,  angle  of  scattering,  or  emerging  energy),  it  is  necessary 
to  know  which  perturbations  are  in  effect  for  the  current  neutron  energy.  For 
this  purpose,  the  entire  energy  range  of  the  problem  is  broken  up  into  energy 
bins  whose  boundaries  are  the  set  of  the  upper  (EH1)  and  lower  (EL1)  limits 
of  all  the  perturbations  i (cf„  Section  3,2).  These  (EH1,  EL1),  i=l,n,  to- 
gether with  the  overall  energy  limits  of  the  problem  EHIGH  and  ELOW  are  then 
sorted  in  decreasing  magnitude  and  entered  into  a table  (the  PETAB-array) . 
Another  table  (IPBIN-array)  is  constructed,  indicating  the  numbers  of  perturba- 
tions in  the  various  bins  of  PETAB-array  and  also  the  ID's  of  these  perturba- 
tions (the  ID's  are  called  IP;  cf.  Section  3.2).  Both  PETAB-  and  IPBIN-arrays 
are  constructed  in  subroutine  PPROCS,  which  is  called  by  subroutine  INPUTP  of 
program  SAMIN„ 

During  later  Monte  Carlo  calculations,  when  a particle  enters  into  colli- 
sion at  energy  E,  the  energy  table  (PETAB)  is  consulted  to  locate  the  energy  bin 
containing  E„  Knowing  this  bin-number,  one  can  obtain  all  the  information  on 
the  perturbations  affecting  energy  E from  the  data  in  the  IPBIN-array. 


IB 


3.  PROGRAM  SAMIN:  NEUTRON  PERTURBATION  INPUT  PROCESSOR 


3.1  General  Discussion 

SAMIN  is  a driver  program  which  calls  subroutine  INPUTP  to  read  in  and 
process  the  neutron  perturbation  input  data  and  which  writes  the  following 
perturbation  information  on  a tape  or  disk  file: 

MISTER-SISTER  array  processed  perturbation  input  data  (refer  to 

storage  map  of  MISTER-array  in  Appendix  C) . 

NTYP-array  cumulative  array  Of  perturbations  by  type,  i.e., 

NTYP (l)=number  of  perturbations  type  1,  NTYP(2)= 

NTYP ( 1 ) tnumber  of  perturbations  of  type  2,  etc. , 
finally  NTYP(N)=sum  of  all  perturbations  of  types 
1,  2,  . . , and  N. 

LOCLIP  location  where  perturbation  data  for  a specific 

element  of  a specific  type  begin  in  MISTER-array. 
(Refer  to  storage  map  of  MISTER-array  in  Appendix  C. ) 

NPROB  total  number  of  problems;  i.e  , unperturbed  plu? 

perturbed  problems. 

LLNEXT  location  of  first  available  storage  +1  after  ele- 

mental perturbation  data  in  MISTER-array.  Sampling 
CHI-  and  ENN-tables  start  at  location  (LLNEXT- 1 ) if 
such  perturbations  exist. 

MMAX  dimension  of  MISTER-array. 

EIDW  low  energy  cut-off  limit  of  problems  run. 

EHIGH  high  energy  limit  of  problems  run. 
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PETAB-array  perturbation  energy  table  where  PETAB (1)=EHIGH, 

PETAB(last)=ELOW,  intermediate  entries  including  all 
the  perturbation  energy  limits  (All  EH's  and  EL's, 
cf.  Section  3.2.) 

IPBIN-array  perturbation  binning  table  (refer  to  storage  map 

description  of  IPBIN-array  in  Appendix  D) . 

LP  total  number  of  energy  bins  in  PETAB-array 

3.2  Processor  INPUTP  and  Description  of  Input  Data  for 

Program  SAMIN 

A Description  of  INPUTP,  the  Routine  for  Reading  and  Processing 

Perturbation  Data  and  Perturbation-Problem  Correspondences 

As  mentioned  previously,  10  types  of  perturbations  are  permitted.  All  of 
the  numerical  data  for  the  perturbations  are  read  in  by  subroutine  INPUTP,  with 
the  exception  of  those  for  perturbations  of  type  2 (whole  element).  For  this 
type,  only  the  element  (nuclide)  identification  integers  for  the  basic  and  per- 
turbing data  are  specified,  and  it  is  assumed  that  both  of  these  identification 
integers  have  been  included  in  the  element  input  to  the  previously  run  cross- 
section-processing program  SAM-X. 

For  many  of  the  input  items  listed  below,  e.g.,  those  numbered  1;  4.1b; 

4. 3-4. 5a, b;  4.6a,b;  etc.,  certain  input  options  have  been  implemented  to  faci- 
litate the  entry  of  single  floating  numbers  or  arrays  of  floating  numbers  with- 
out repeating  the  right- justified  E+nn  punch  on  the  input  cards. 

Specifically,  the  first  of  these  options  permits  the  cards  calling  for 
the  EELOW  and  EEHIGH  of  the  problem  (card  1)  and  all  the  perturbation  cards 
calling  for  an  EL  and  EH  (e  g.,  card  4.3a)  to  have  the  (e.g  ) E+06  punch  omitted, 
and  replaced  by  *6  in  columns  79  and  80  of  the  pertinent  cards.  The  code  then 
enters  EL  and  EH  in  memory  as  the  actual  constants  of  the  designated  data  fields. 
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multiplied  by  10  . In  general,  for  such  cards  (and  others  indicated  below) 

+N 

the  pertinent  input  quantities  are  scaled  (prior  to  storage)  by  10  if  *N 

-N 

appears  in  columns  79  and  80,  and  by  10  if  /N  appears  in  these  columns. 

This  scaling  is  done  by  subroutine  SCALE. 

Cards  for  which  this  option  is  provided  are  indicated  in  the  following 
data  and  format  description  by  the  entry  *N  at  the  right-hand  side  of  the 
format  statement,  the  * is  to  be  understood  as  representing  * or  /,  while  N 
an  integer  or  blank;  (a  non-integer  alphanumeric  in  column  80  will  cause  a 
system  error  since  the  code  scans  this  column  in  II  format).  If  column  80  has 
a non-zero  numerical  entry,  the  code  will  do  nothing  to  the  quantity  in  the 
scaled  field  unless  column  79  contains  an  * or  a /.  In  the  input-data  format 
listings  below,  the  quantities  affected  by  the  optional  scaling  are  underlined. 

Two  other  types  of  optional  scaling  of  input  data  are  provided 

The  first  of  these  applies  to  single  arrays  of  floating  numbers  such  as 
energies  (as  in  items  4.6b,  4.7c).  For  most  of  these  data,  called  for  in 
7E11  4 format,  *N  or  /M  in  columns  79  and  80  of  the  first  card  of  the  array  will 
cause  all  the  input  numbers  to  be  multiplied  (prior  to  printing  and  storing) 
by  10+N  or  10-M  respectively;  an  optional  P (or  any  other  non-blank  character) 
in  column  78  will  cause  the  data  in  the  array  to  be  printed  out  in  a highe’ — 
precision  format  (6E16.6)  instead  of  the  normal  or  default  format  (10E12.4). 

This  scaling  and  printing  is  done  in  subroutine  READ75. 

Arrays  for  which  this  option  is  provided  are  indicated  below  by  the  entry 
P*N  at  the  right-hand  side  of  the  page;  P stands  for  any  non-blank  character, 

* stands  for  * or  /,  and  N is  the  integer  for  exponent  scaling. 

Finally,  there  are  several  types  of  data  which  are  read  in  as  coordinated 
tables  (6E11  4,  x and  y alternating;  as  in  energy,  cross  section,  items  4 3- 
4.5b,  4 7g ; or  energy,  probability,  as  in  item  4 7i).  For  most  of  the  input 
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data  of  this  type,  columns  77-80  of  the  first  card  in  the  pertinent  read 


1 


operation  may  have  entered  a notation  such  as  *L*N,  where  each  * may  indepen- 
dently be  * or  /,  and  L and  N are  integers „ Entries  in  columns  77,  78  scale 
the  x (first,  or  independent-variable  field)  and  those  in  columns  79,  80 
affect  the  y (second,  or  dependent- variable  field).  Hence  E,  cards  punched 
as  follows: 

Col.  Col.  Col.  Col.  Col.  Col 

11  22  33  44  55  66  77  78  79  80 

Card  1 1.0  7.1  1.5  6.1  2.  5.1  ...»  6 / 3 

Card  2 2.5  4.1  3.0  3.1  (blank) 

could  cause  energies  1.0E+6ev,  ..,  3.0E+6ev  to  be  correlated  with  cross  sec- 
tions of  .0071,  .0061,  . . , 0031  barns,  respectively,  in  portions  of  the  INPUTP 

code  where  such  scaling  is  provided. 

Either  or  both  of  columns  78,  80  may  be  blank,  but  neither  should  contain  a 
non-numeric  punch.  Scaling  action  is  taken  on  a field  only  if  the  operator  (in 
col.  77  for  x or  79  for  y)  is  * or  /.  This  scaling  is  done  in  subroutine  RCTAB. 
Quantities  which  may  be  scaled  by  using  this  option  are  underlined  in  the  format 
lists  below. 

INPUT  to  INPUTP 

The  following  describes  the  input.  The  information  is  summarized  in 

Section  3.3 

Item  1 - EE LOW,  EEHIGH,  ISSW  (1,  2,  ...,  12);  FORMAT  (2E11.4,8X,12I1) , *N 

The  first  two  quantities  are  the  low-energy  cut-off  and  the  high-energy 

* 

limit,  both  in  ev,  for  tracking  in  the  main  Monte  Carlo  calculation  . 

ISSW  (1,2, ,12)  are  "sense  switch"  digits  used  for  selective  debugging 

and  additional  information  via  printouts.  Non-zero  entries  elicit  the  printouts. 
At  present  only  indices  1,2,6,  and  11  have  any  effect.  These  are: 

*See  Appendix  H 
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ISSW(l):  Prints  out  computed  Legendre  moments  of  tabulated  angular 


distributions . 

ISSW(2):  Prints  out  entire  perturbation-array  MASTER*  at  end  of 

calculation,  in  several  formats. 

ISSW(6):  Prints  out  details  of  ANTERP  routine  operation,  locating 

abscissa  corresponding  to  specified  ordinate  (area)  in  elastic 
angular  distribution. 

ISSW(ll):  Prints  out  partial  MASTER*  array,  showing  unedited 

problem-to-perturbation-correspondence  data.  These  data  are  also 
edited  and  printed  independently  of  ISSW(ll). 

Item  2 - NTYP  (1 , 2 , . . . , 10) ; FORMAT  (1015). 

Ten  integers,  indicating  the  numbers  of  perturbations  present  of  types 

1,2,..  ,10,  respectively.**  Any  or  all  of  these  may  be  zero.  If  all  ire  zero, 

this  card  is  the  last  one  that  is  required.  Otherwise: 

10 

Item  3 - For  each  perturbation  KK,  where  KK=1,2,...,  j NTYP(J),  a card  wit' 

J=1 

KK,  NPROB,  (JPROB(L)  ,L-1, NPROB)  ; FORMAT  (1215) 

Here  KK=perturbation  number;  NPROB  is  the  number  of  output  problems  to  be 
affected  by  the  perturbation  (CKNPROB<10) ; and  JPROB (L)  are  the  particular 
problem  numbers  (all  distinct,  and  each  JPROB^IO) . KK  must  be  present;  if 
NPROB=0,  subsequent  input  for  that  perturbation  is  limited  to  a single  card. 

* NOTE:  Not  to  be  confused  with  MASTER  array  in  Program  SAMCAR. 

**  The  array  is  internally  transformed  to  a cumulative  one,  to  conform 
to  the  definition  of  Section  3 1. 
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These  data  are  followed  by  cards  describing  the  individual  perturbations, 


as  follows. 

All  data  for  perturbations  of  type  1 (if  any)  must  be  given  first,  and 
numbered  (KK)  in  order  from  1 to  NTYP(l);  those  for  type  2 (if  any)  must  follow, 
and  be  numbered  consecutively  from  NTYP(1)+1  to  (NTYP(l)  + NTYP(2));  etc„ 

Header  cards  must  be  present  for  each  perturbation,  with  identifiers  KK=1,2,..., 

10 

Z NTYP(J),  including  those  for  which  NPROB=0  (item  3,  above). 

J=1 

For  each  KK  with  NPROB=0,  no  further  data  are  required.  Otherwise,  the  required 
inputs  are  as  follows:  (Items  4.1  - 4.10). 

Item  4.1  - Type  1,  Composition  (Present  if  NTYP(1)>0). 

Here  a specific  composition  (set  of  concentrations),  one  of  those  processed 
by  BAND,  is  perturbed.  Required  input: 

a)  KK,  ICOMP,  NLM;  FORMAT  (X, 14,215), 

KK=perturbation  index;  (an  asterisk  or  other  mark  may  be 
punched  in  column  1 to  increase  legibility  of  the  data  deck); 

ICOMP  = composition  number,  as  specified  for  BAND; 

NLM  = number  of  elements  in  composition,  as  specified 
for  BAND;  followed  by 

b)  CONC (1,2,, ■ . , NLM) ; FORMAT  (7E11.4),  P*N. 

These  are  the  NLM  concentrations  of  the  elements  comprising 
the  composition,  in  atoms/ (barn-cm) , given  in  the  same  order 
as  provided  in  the  BAND  input. 

NOTE:  In  all  the  following  sections  (4.2-4.10),  it  is  to  be  understood 

that  for  any  perturbation  KK  for  which  NPROB=0,  no  data  past  the  first  card  are 
to  be  supplied. 
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Item  4,2  - Type  2,  Whole  Element  (Present  if  NTYP(2)>0) 


Here,  a specified  base-case  element  is  understood  to  have  all  its  cross 

section  data  perturbed  in  the  entire  energy  range.  The  perturbed  data  are  thos^ 

★ 

of  another  specified  element  present  on  the  existing  element  data  tape  (ED T) . 
Input: 

a)  KK,  ID1,  ID2;  FORMAT  (X, 14, 2110). 

KK  = perturbation  index; 

ID1,  ID2  are  base-case  and  perturbation  element  integer  identifiers, 
e.g„,  92235  and  92000.  Conventionally  an  identifier  is  of  the  form 
IZAAA,  where  IZ  is  an  atomic  number  and  AAA  an  atomic  weight.  AAA=000 
is  frequently  used  to  indicate  a natural  element  (mixture  of  isotopes) 

Items  4.3  - 4.5  - Types  3,4,5:  Perturbation,  within  a specified  energy 

range,  of  total  (ct),  scattering  (o  ),  or  inelastic  (a.  ) 

t s in 

microscopic  cross  sections,  respectively,  of  a specified  element. 

Here,  it  is  implied  that  specifying  a perturbation  of  a alone  implies  a 
change  in  the  absorption  cross  section;  explicitly  changing  og  alone  means 
changing  the  elastic  and  absorption  cross  sections;  and  explicitly  changing 
cr  alone  means  changing  the  inelastic  continuum  and  the  elastic  scattering 
cross  sections.  These  relationships  can  be  understood  from  Figure  1,  bei  ,w, 
which  indicates  the  components  of  the  total  microscopic  cross  section;  the  com- 
ponents in  boxes  are  those  which  are  specified  in  the  basic  or  perturbation  data 


As  the  BAND  tape  is  prepared  independently  of  the  present  input  processing, 
it  is  required  that  the  BAND  input  includes  the  perturbing  element.  This 
can  be  achieved  by  declaring  all  perturbing  elements  to  be  constituents  of 
a dummy  composition  (i„e.,  one  which  is  not  referenced  by  any  physical  region) 


(The  level-excitation  cross-section 
o,  is  actually  supplied  as  cross- 
sections  for  the  excitation  of  the 
individual  levels.) 


Figure  1 - Components  of  Microscopic  Cross-Section 


By  a suitable  superposition  of  perturbations  of  types  3,4,  and  5 it  is  possible 
to  specify  any  desired  perturbation.  Thus,  for  example,  if  it  were  desired  to 
increase  a while  leaving  o unchanged,  this  would  require  simultaneously  speci- 

G cl 

fying  the  required  changes  (as  functions  of  neutron  energy)  in  (type  4)  and 
a (type  3). 

Required  input  (types  3,  4 and  5): 

For  each  perturbation  present  of  these  types: 

a)  KK,  NEP,  JZ , EL,  EH,  FMULT ; FORMAT  (X, 14, 15, 110, 3E15. 5) , *N 
KK  = perturbation  index; 

NEP  = number  of  energy  points  at  which  the  cross  section  is  specified 
(at  least  2 for  tabulated  data; for  NEP<2,PMULT  is  used  as  the 
multiplier  to  be  applied  to  the  unperturbed  cross  sections  in  the 
energy  range  from  EL  to  EH  to  obtain  the  perturbed  values); 

JZ  = element  identifer  as  given  by  SAM-X  processor; 

EL,  EH  = low  and  high  energy  limits  (ev)  over  which  the  perturbation 


is  to  apply; 
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FMULT  = multiplier  for  cross  sections  (used  if  NEP<2); 


next,  if  NEP>2  (omit  (b)  if  NEP<2) . 

b)  Supply  NEP  pairs  of  (energy  (ev) , cross-section  (barns)): 

E (L) , 0(L) , for  L=1 , 2 , . . . , NEP ; FORMAT  (6E11 . 4 ) , *L, *N 

The  energies  and  corresponding  cross  sections  are  given  alternately, 

3 pairs  to  a card,  energies  in  increasing  order.  As  in  all  the  per- 
turbations involving  a specified  energy  range,  any  supplied  table  of 
energies  must  span  the  perturbation's  EL  and  EH;  i.e.,  E(l)  <EL  and 
E (NEP)  >EH,  otherwise  the  problem  stops.  The  spacing  of  the  energy 
points  should  be  such  that  linear  interpolation  of  a in  E is  permissible 

Item  4 6 - Type  6,  Angular  Distribution  of  Elastic  Scattering 
Required  Input: 

a)  KK,NE,ID,EL,EH;  FORMAT  (X , 14 , 15, 110 , 2E1 5. 5 ) , *N 
KK  = perturbation  index; 

NE  = number  of  energy  points  at  which  perturbed  distributions 
are  supplied; 

ID  = element  identifier; 

EL,  EH  = low  and  high  energy  limits  (ev)  of  the  perturbation 

b)  E (L) , L=l,  NE  (energies  (ev)  in  increasing  order,  for  which 
distributions  are  given);  FORMAT  (7E11  4),  P*N,  (E(l)d:L,  E (NE) >EH) . 

The  spacing  of  the  energies  should  be  such  that  linear  interpolation 
in  energy  of  the  implied  distributions  is  permissible.  Next  (c-c), 
CHI-table  criteria;  one  package,  applicable  for  all  of  the  NE  energies: 
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c)  LADHOC,  N;  FORMAT  (215) 


LADHOC,  Maximum  allowable  number  of  CHI's  (equi-probable 


1-U 

cm 

values  of  \ = — , exclusive  of  the  implied  x=0  and  1), 

and  N,  number  of  Legendre  moments  of  the  computed  X-table  to  be 
compared  with  those  of  the  supplied  distribution. 

LADHOC<40,  N<_50„ 

If  LADHOC  = 0;  default  values  will  be  set  internally  for 
LADHOC  (=30),  N (=3) , e (=.02,  .02,  .02)  and  a(=.5); 

X— 1 , 2 , j 


where  the  e.  and  a are  defined  below, 
x 

Skip  inputs  (d,e) , and  go  to  f. 

If  LADHOC>i:  Let  N be  the  value  supplied  on  card  c;  then  if  N=0, 
skip  to  (e) (default  values  of  N and  £ will  be  used) . 

If  N^O,  supply  (d)„ 

d)  £(I),  1=1, N FORMAT  (10F8.3)  (no  scaling  allowed) 

Blank  entries  for  e,  or  values  <„001,  are  each  replaced  internally 
by  001..  The  e are  absolute  deviations  to  be  allowed  in  the  com- 
puted Legendre  moments  of  the  x-table „ 


Moment 


f(y)P.  (U)dvi,  where  f(p)=^of  the  supplied  (and 
■>  dli 


3 


code-normalized)  angular  distribution,  and  IV  00  is  the  jth 
Legendre  polynomial.  The  implied  moments  refer  to  those  beyond 
the  zero-th,  which  is  unity. 


e)  a FORMAT  (E12.4) 

a is  the  fractional  deviation  to  be  allowed  in  the  dp/dy  implied 
by  the  generated  x-table,  as  compared  with  the  maximum  value  of 
the  supplied  distribution,  within  any  one  X”bin.  (Actually,  the 
value  of  the  distribution  within  the  x-table  must  deviate  from 
the  bin-maximum  supplied  value,  and  from  the  average  value  of  the 
distribution  (=  5),  by  more  than  a in  order  for  a particular 
X-table  to  be  rejected.  The  e-tests  are  the  final  criteria  to  be 
applied.)  The  code  tries  x~tables  with  more  and  more  CHI's  (up 
to  LADHOC)  to  attempt  to  satisfy  all  criteria. 

f)  LTT;  FORMAT  (12) 

Type  of  data  supplied,  LTT=1  (Legendre  coefficient  expansion)  or 
LTT=2  (tabulated  y, dp/dy). 

The  choice  of  LTT  fixes  the  type  of  distribution  input  for  all 
the  input  energies. 

Next;  If  LTT=2,  skip  to  Item  4.6„2. 

Item  4.6.1  - applies  if  LTT=1,  supply  items  (g),  (h)  for  each  of  the 
NE  energies,  then  skip  to  Item  4 7. 

(g)  NL,  Number  of  coefficients  (past  the  zero-th)  to  be  used  in  the 
Legendre  expansion  of  the  normalized  center-of-mass  angular  distri- 
bution f (y ) . 

(NL<50)  FORMAT  (15) 

(h)  c (L) , L=l,  NL;  FORMAT  (7E11  4)  (no  scaling  allowed) 
values  of  expansion  coefficients, 

-/+1 

c 1/  I f(y)Pf(y)dy  of  the  normalized  f(y),  where 
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2k+l 

'<">  ■ ~r  W”  v1-0 

k=0 

Pq(m)  = 1;P1(U)  = V ; P 2<U)  = ~ ; 

mP  (y)  = (2m-l)uP  AV)  - (m-l)P  _(W). 
m m- 1 m-  2 

Item  4.6.2  - applies  if  LTT=2 . Supply  items  (i),  (j)  for  each  of  the 
NE  energies,  then  enter  Item  4.7. 

(i)  MNP,  number  of  (h, dp/dP)  pairs  over  range  -l^U^+1  supplied 
to  describe  center-of-mass  angular  distribution  for  particular 
energy.  (2^MNPflOO) ; FORMAT  (15) 

j)  (y,dp/dy),  (MNP  pairs);  FORMAT  (6E11.4)  (no  scaling  allowed) 
center-of-mass  cosine  of  scattering  angle  and  probability  per 
unit  cosine  (differential  cross  section  or  probability,  not 

necessarily  normalized);  y and  dp/dy,  alternately,  3 pairs  to  a 

I . 

card  The  first  dp/dy  should  correspond  to  y=-lo,  the  last  to 
y=+l . ; (the  code  sets  U ( 1 ) =- 1 . and  y(MNP)=+l.)„  The  spacing 
should  be  such  that  linear  interpolation  of  dp/dy  in  y is  permiss- 
ible. (Repeat  (i)  and  (j)  for  each  subsequent  energy.) 

End  of  Type  6 input. 
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Item  4„7  - Type  7,  Secondary-Energy  Distribution  of  Inelastic 
Continuum  Scattering 

Input : 

a)  KK,NES, ID, EL, EH;  FORMAT  (X,I4,I5,I10, 2E15.5) , *N. 

KK  = perturbation  index; 

NES  = number  of  energy  points  at  which  perturbed  distributions 
are  supplied; 

ID  = element  identifier; 

EL, EH  = low  and  high  energy  limits  (ev)  of  the  perturbation 

b)  NENN(L),  L=1  NES;  FORMAT  (1615) 

These  integers  are  the  numbers  of  ENN  (secondary-energy 
boundaries,  defining  equal ly-probably  emergent-energy  bins)  to  be 
generated  for  each  of  the  NES  incoming  energies.  The  values  of  NENN 
should  be  >2,  and  non-decreasing  (the  code  corrects  the  list  if 
necessary  to  make  NENN(1)>2  and  d (NENN)/dL_>0)  . 

c)  E (L) , L=1 , NES;  FORMAT  (7E11.4) ,P*N. 

NES  values  of  incoming  energy,  low  to  high.  The  list  should  at 
least  span  (EL, EH). 

d)  NREAC,  Number  of  reactions  contributing  to  the  inelastic 
continuum.  FORMAT  (15) . 

Supply  Package  (e)-(n)  (with  appropriate  internal  choices)  for 
each  of  the  NREAC  reactions. 

e)  KREAC,  MT,  NK;  FORMAT  (X,I4,2I5). 

Here  KREAC  is  the  (required)  reaction  number  (1,2, . . . , NREAC) ; 

MT  is  the  ENDF/B  identifier  for  the  reaction  type;  the  allowed 
values  of  MT  and  their  meanings  are: 


J1 


MT 

Meaning 

4 

inelastic  (n,n') 

16 

(n, 2n) 

17 

(n, 3n) 

22 

(n,nh  ) 

23 

(n,n'3a) 

24 

(n,  2na) 

25 

(n,  3na) 

28 

(n,n'p) 

(use  of  any  other  MT  will  stop  the  problem);  NK  is  the  number  of 
supplied  energy  distributions,  contributing  to  the  total  energy  dis- 
tribution for  reaction  KREAC. 

f)  NEPS,  number  of  energy  points  at  which  the  cross  section 
for  reaction  KREAC  will  be  supplied. 

(NEPS>2);  FORMAT  15) 

g)  (ES (L) , 0k(L),  L=l,  NEPS);  FORMAT  (6E11 . 4 ) , *L*N 

NEPS  pairs  of  (energy,  cross  section  for  this  reaction);  three 
pairs  to  a card. 

(The  energies  should  be  increasing,  and  spaced  so  that  linear  in- 
terpolation for  is  permissible;  the  span  of  energies  for 
need  not  encompass  EL,  EH;  for  an  energy  E(L)  which  lies  outside 
the  range  of  energies  for  which  is  defined,  the  value  of  is 
taken  to  be  zero.) 

Within  this  reaction-package,  supply  (h-n)  for  each  of  the  NK  distributions; 

h)  NEPP,  Number  of  energies  at  which  the  probability  for  this  distri- 
bution (fraction,  p,  of  to  be  associated  with  this  distribution)  will 
be  specified  (NEPF>2);  FORMAT  (IS). 


i)  (EP  (L)  , P (L)  , L=1 , NEPP);  FORMAT  (6EU.4)  ,*L*N 


Alternate  energy,  probability  for  this  distribution;  3 pairs  to 
a card,  energies  increasing.  Outside  the  tabulated  range  of  EP, 

P is  assumed  to  be  zero. 

Note:  For  any  one  reaction  at  a specified  energy  E,  the  sum  of 

the  p's  for  the  various  distributions,  as  obtained  by  interpolations 
in  E,  should  equal  1: 


NK 

E p. (E)  = 1. 
i=l 


The  code  does  not  check  this,  so  the  problem  originator  should. 

(Allowing  Ep.  to  be  different  from  1 means  that  the  generated 
i 

energy  spectrum  is  not  weighted  correctly  by  o^,  but  effectively  by 

a,  » Ep . ) o 
k *i 

The  satisfying  of  this  condition  (EP£(E)=1  ) is  facilitated  if  the 

i 

same  energy  mesh  is  given  for  all  of  the  tabulated  p(E)„  Then  if 
the  condition  is  satisfied  by  the  p's  at  each  energy  nesh  point,  it  will 
be  satisfied  by  the  linearly-interpolated  p's  at  al i intermediate  energies, 
j)  LF,U+; FORMAT  (I 5, E15 „ 4) , *N 

LF  is  the  energy  distribution  type;  LF=5  or  9;  U ( >0  ) is  pertinent  if 

LF=9;  it  is  used  to  truncate  an  otherwise  Maxwellian  energy  distribution 

at  a maximum  E'  of  E-U.  The  meanings  of  the  LF's  are  as  follows: 

LF=5:  Spectrum  is  supplied  as  a single  table  of  X,g(X),  where  X=E'/0, 

0 is  0(E),  a nuclear  temperature,  specified  as  a tabulated  function  of  E, 

and  g (X) =dn/dX'Vj2.,  the  probability  density  of  the  secondary-energy  spectrum; 
qE 

LF=9:  Spectrum  is  Maxwellian, 


e M,  K Drake  (ed.),  "Data  Formats  and  Procedures  for  the  ENDF  Neutron  Cross 
ction  Library",  BNL  50274,  ENDF  102,  Vol„  1 (1970). 


where  U is  constant  for  all  incoming  energies  E,  and  0 is  a tabu- 
lated function  of  E. 

The  code  generates  the  required  normalization  for  both  cases  LF=5 
and  LF=9,  so  that 


If  a particular  table  of  g(X)  (LF=5)  does  not  extend  to  high  enough 
X for  a particular  E to  include  E'=E,  the  code  sets  g(X)=0  over 
the  undefined  range „ 

k)  NEPTH,  Number  of  energy  points  at  which  9(E)  is  supplied., 

(NEPTH>2);  FORMAT  (15) 

l)  (ET (L) , 9 (L) , L=l, NEPTH);  FORMAT  (6E11 . 4) , *L*N 

Alternate  energy,  theta  for  this  distribution.  3 pairs  to  a card. 

Energies  increasing,  spaced  to  permit  linear  interpolation  for  theta. 

For  an  incident  E outside  the  9-defining  energy  range,  g(X)  or  dn/dE' 
is  set  equal  to  zero  for  all  E',  O^E'^E. 

m)  (Supply  only  if  LF=5;  skip  (m)  and  (n)  if  LF=9) ; 

NPX,  Number  of  values  of  X(=E'/0)  at  which  g(X)  is  supplied 
(NPX>2);  FORMAT  (15) 

n)  (X (L) , g (L) , L=l,  NPX)  ; FORMAT  (6E11 . 4) , *L*N 

Alternate  X,  g(X)  for  this  distribution.  3 pairs  to  a card  In- 
creasing X.  At  any  E'  with  X=E'/® (E)  outside  the  tabulated  range, 
g(X)  is  taken  to  be  zero.  Linear  interpolation  of  g(X)  in  X is 
assumed. 

Repeat  items  h-1  (LF=9)  or  h-n  (LF=5)  for  each  subsequent  distribu- 
tion 2,.„.,NK. 

Repeat  items  e-n  for  each  subsequent  reaction  KREAC=2, . . . ,NREAC. 

End  of  type  7 input. 

J 
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Item  4.8  - Type  8,  Cross  Section  for  Level  Excitation 


Here  the  individual  cross  sections  for  selected  levels  of  a particular 
element  are  specified,  as  functions  of  energy  °n  a single  energy  mesh  for  thi 
perturbation. 

Input: 

a)  KK,NES , ID, EL, EH;  FORMAT  (X, 14 , 15, 110, 2E15. 5) , *N 
Perturbation  index,  Maximum  number  of  energies  at  which  o-levels 
are  given,  element  identification,  perturbation  E-low,  E-high;  : 

b)  NLEVP,  the  number  of  levels  perturbed  (1<NLEVP<L00) ; 

FORMAT  (15) 

C)  (KLEV(K),  NEL(K),  K=l,  NLEVP);  FORMAT  (1615) 

KLEV (K)  is  the  level  number  of  the  Kth  level  perturbed; 

NEL(K)  is  the  number  of  consecutive  energies  in  the  input  descend- 
ing-energy mesh  (item  (d) , below)  at  which  the  cross  section  for  level 
KLEV  is  supplied  (from  the  first  (highest)  energy  downward) . Alter- 
nate KLEV,  NEL,  8 pairs  to  a card. 

NOTE:  1.  The  KLEV's  must  be  increasing,  but  need  not  be 

consecutive  integers, 

2 The  first  NEL  must  be  ^NES, 

3.  The  NEL's  must  be  non-increasing,  and 

4.  The  minimum  value  of  NEL  must  be  >2„ 

If  these  conditions  are  not  satisfied  by  the  data,  the  problem  stops. 

A cross  section  must  be  defined  for  at  least  two  energies  since  outside 
the  range  of  energie:.  for  a cross  section  the  perturbation  is  taken 
to  be  non-existent,  i.e.,  interpolation  is  not  done  between  a perturbe  ; 
and  an  unperturbed  cross  section. 
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d)  E (L) , L=1 , NES;  FORMAT  (7E11.4),P*N 


Item 


energies  at  which  level  cross  sections  are  specified;  in  decreasing 
order;  E(1)_>EH  for  this  perturbation,  and  E(NES)_<EL. 

The  energy  mesh  should  be  chosen  so  that  for  any  perturbation  level, 
the  cross  section  can  be  supplied  for  at  least  2 energies,  and  such 
that  linear  interpolation  in  energy  is  correct, 
e)  For  each  perturbed  level,  supply 
(^(.1)^=1,  NEL(K);  FORMAT  (7E11.4);P*N 

Cross  Sections  corresponding  respectively  to  the  first  (highest)  NEL(K) 
energies  of  the  list  (d).  7 to  a card. 

Start  a new  card  for  each  new  level.  If  scaling  is  used,  one  must  enter 
P*N  for  the  first  card  of  each  new  level.  Scaling  may  differ  from 
level  to  level. 

4.9  - Type  9,  Angular  Distribution  of  Discrete- Level 
Inelastic  Scattering 

The  first  portion  of  the  data  required  (items  (a)-(d))  is  the  same  as  for 
type-8  input;  restrictions  on  the  NEL(K)  are  less  stringent. 

Input: 

a)  KK,NES,ID,EL,EH;  FORMAT  (X, 14 , 15, 110, 2E15. 5) , *N 

NES  = number  of  energies  at  which  (y,dp/dy)  tables  are  given  for 
selected  levels;  NES^2. 

b)  NLEVP,  number  of  levels  perturbed  (l£NLEVP£l00) ; 

FORMAT  (15) 

c)  (KLEV(K),  NEL(K),  K= 1 , NLEVP);  FORMAT  (1615) 

KLEV(K)  is  the  level  number  of  the  Kth  level  perturbed, 

NEL(K)  is  the  number  of  consecutive  energies  in  the  input  descending- 
energy  mesh  (item  (d) ) at  which  angular-distribution  tables  (p,dp/d;i) 
for  level  KLEV  is  supplied,  from  the  first  (highest)  energy  downward. 

8 pairs  to  a card. 


.36 


NOTE: 


1.  the  KLEVS  must  be  increasing,  but  need  not  be 


consecutive  integers, 

2 2^NEL(K)^NES,  and 

3.  the  NEL's  are  not  otherwise  constrained. 


d)  E (L) , L=l,  NES;  FORMAT  (7E11.4),  P*N 

energies  at  which  angular-distribution  tables  are  specified;  in 
decreasing  order;  E(1)^EH,  and  E(NES)^EL„  The  energies  should  be 
close  enough  to  permit  linear  energy  interpolation  for  dp/dy  at  a 
given  y. 

For  each  perturbed  level,  supply  items  (e): 

e)  Within  this  level,  supply  e.l,  e.2  for  NEL  (level)  energies. 

(e.l)  NMUP ( >2 ) ; FORMAT  (15) 

No  of  (y,dp/dy)= (cos9  ,dp/dy)  pairs  supplied  over  range  -l<y<+l 

c .in. 

(e.2)  y (L) , dp/dy  (L) , L=l,  NMUP;  FORMAT  (6E11 . 4 ) , *L*N 
Three  pairs  to  a card.  The  first  dp/dy  should  correspond  to  y=-l,  the 
last  to  y=+l;  (the  code  sets  y(l)=-l  and  y(NMUP)=+l)  The  values 
of  dp/dy  should  correctly  represent  the  shape  of  the  desired  distri- 
bution but  they  need  not  be  normalized;  this  is  done  by  the  code 


so  that 


1.  Linear  interpolability  of  dp/dy  in  y is 


assumed. 

(Supply  e.l,  e.2  for  all  energies  E,  (highest)  to  E.^T  for  this  level.) 

1 NEL 

(Repeat  package  (e.l,  e.2)  x NEL^eve^  f°r  each  successive  perturbed 
level . ) 

End  of  type-9  input 
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Item  4„10  - Type  10,  Source  Spectrum 
Input: 

a)  KK, NES;  FORMAT  (X,I4,I5) 

Perto  index;  No„  of  energy-boundaries  specified  for  source-spectrum 
histogram  (NES>2) „ 

b)  (E (L) , S (L) , L=1 , NES);  FORMAT  (6Ell„4) , *L*N 
NES  pairs  of  (Energy,  As/AE);  energies  increasing; 

E ( 1 ) <EEIX3W  of  entire  problem  (first  card  of  input),  and 

AS 

E (NES) >EEHIGHC  The  source  density  S(L)=  — represents  the  constant 

AE 

value  referring  to  the  range  from  E(L)  to  E(L+1);  hence  S(NES-l)  is 
the  last  pertinent  value  of  S read„ 

This  concludes  the  perturbation-data  inputs  The  following  section  gives  a 
resume  of  the  input  and  formats,. 
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3.3  Summary  of  Input  Data  for  Program  SAMIN 


Item 


Item 


Item 


Item 


Item 


DATA  and  Restrictions 
1„  EELOW,  EEHIGH,  ISSW (1 , 2, . . 0 , 12 ) 

Energy  limits  for  tracking,  and 
sense  switch  options 
2.  NTYP(1,2,... ,10) 

NTYP (K)=number  of  perturbations  of 
type  K present  in  the  data.  (If 
all  are  zero,  end  of  input). 

Otherwise: 

30  For  each  perturbation: 

KK,NPROB,  (JPROB(L),  L=1,NPR0B) 

KK=pert0  noo,  NPROB=no.  of  affected 
problems  (<10) ; 

JPROB (1, 2, o . . ,NPROB)  are  distinct 
numbers,  each  <10.  NPROB  may  be  zero; 
if  it  is,  data  for  pert.  KK  is  subse- 
quently limited  to  1 card  (that  with 
item  KK  in  columns  2-5) „ 

10 

Data  for  specific  perturbations  with  KK=1,2,0q.,  I NTYP (L) : 

L=1 

401  - Type-1,  Composition 

a)  KK,  ICOMP,  NLM  (X, 14, 215) 

Pert,,  composition  no0,no„  of  elements 

b)  CONC  (1,2,..., NLM)  (7E11„4),P*N 

Concentrations  (atoms/barn-cm) 

4 2 - Type  2,  Whole  Element 

a)  KK,  ID1,  ID2  (X, 14, 2110) 

Pert.,,  base  ID,  perturbing  ID„ 


FORMATS 

(2E11.4,8X,12I1) ,*N 
(1015) 


(1215) 
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(XrI4, 15,110, 3E15. 5) , *N 


Items  4,3-4  5 - Types  3,4,5  (o  ,o  a.  ) 

' ""  t S in 

a)  KK,  NEP,  ID,  EL,  Etl,  FMULT 
Pert.,  no.  of  energy  pts.,  element 
ID,  E-low,  E-high,  multiplier  for 
X-sects  (used  if  NEP<2). 

b)  (omit  b if  NEP<2) 

If  NEP-^2 : 

(E,o)  for  L=1 , 2 , . . . , NEP  (6E11 . 4 ) , *L*N 

(Energy  (eV) , cross  section 
(barns))  pairs,  three  to  a card. 

Energies  increasing,  must  span  EL 
and  EH. 

Item  4.6  - Type  6,  Angular  Distribution  of  Elastic  Scattering 

a)  KK,  NE,  ID,  EL,  EH  (X , 14 , 1 5 , 110 , 2E15 . 5 ) , 

Pert  , energies  for  distributions, 
element  ID,  E-low,  E-high 


*N 


(NE>2) 

b)  E ( 1 , 2 , . . . , NE ) (7E11 . 4 ) ,P*N 

Energies  (eV,  increasing); 

must  span  EL,  EH. 

c) -(e):  CHI-table  criteria: 

LADHOC(<40),  N (< 50)  (215) 


[lax.  no.  of  CHI  boundaries, 

No.  of  moments  to  compare 
If  LADHOC=0,  skip  (d)  and  (e); 


(go  to  (f ) ) 

If  LADHOC^l , and  if  N=0,  skip  to  (e) ; 
(if  ^0,  go  to  (d)  ) 


A 


40 


(10F8. 3) 


d)  e(l,2,...,N) 

Allowable  deviations  in  N 
moments  past  the  zero-th 

e)  a 

Allowable  fractional  deviation 
of  x~table  (dp/dp)  for  bin- 
maximum  dp/d(i  tabulated 

f)  LTT  = 1:  Legendre  coefficients 

suppl ied , or 

2:  Tabulated  p,  dp/dp  supplied 
Next  supply  (g) , (h)  if  LTT=1  or  (i),  (j)  if 
LTT=2: 

(LTT=1) : (g)  NL  (£50) 

No.  of  Legendre  coefficients  supplied 
(omitting  zero-th) 

(h)  C (1, 2, . . . ,NL,  Legendre 
coefficients  of  normalized  C.M.  angular 
distr ’ buiion. 

(Repeat  (g),  (h)  for  each  subsequent  energy) 
or 

(LTT=2)  : (i)  MN1  (2jWJP£100) 

No.  of  (p, dp/dp)  pairs  supplied 
j)  (p, dp/dp)  (MNP  pairs) 

First  and  last  dp/dp  should  correspond 
to  if=-l.  and  ii=+l.  respectively;  center- 
of-mass  distribution,  not  necessarily 
normalized. 

(Repeat  (i),  (j)  for  each  subsequent  energy). 


(E12.4) 


(12) 


(15) 


(7E11.4) 


(15) 


(6E11.4) 
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bution  data  are  supplied;  (>2) 

b)  NENN  (1,2, o.o ,NES)  (1615) 

No„  of  ENN  boundaries  for  each  incoming 

energy o Should  be  an  increasing  sequence; 

NENN(1)>2„ 

(Code  corrects  sequence  if  necessary.) 

c)  E (l,2,o«o,NES)  (7Ello4),P*N 

Incoming  energies  (eV)  increasing; 

must  span  EL,  EH 

d)  NREAC  (>1)  (15) 

No.  of  reactions  contributing  to  in- 
elastic continuum.  Supply  items  (e)-(n) 

for  each  reaction  present: 

e)  KKEAC,  MT,  NK  (X, 14, 215) 

Reaction  no. , reaction  type,  no.  of 

distributions  for  this  reaction.  MT  must 
have  one  of  the  values  listed: 


MT 

Meaning 

4 

(n,n' ) 

16 

(n,2n) 

17 

(n,  3n) 

22 

(n,n'a) 

23 

(n,n' 3a) 

24 

(n,2nu) 

25 

(n,  3n«) 

2H 

(n,n'p) 

4 2 


f)  NEPS 


(>2) 


No.  of  energy  points  for  reaction 
cross  section 

g)  (E,ak),  NEPS  pairs  (6E11.4 ) , *L*N 

(Increasing  energies);  assumed  zero 

outside  of  defining  range  of  E. 

Within  this  reaction,  supply  items  (h)-(l) 
or  (h)-(n)  for  each  of  the  NK  distributions. 

h)  NEPP  (>2)  (15) 

No.  of  energies  for  p of  distribution. 

i)  (E,p) , NEPP  pairs  (6E11.4) ,*L*N 

(Increasing  energies) 

j)  LF,U  (I5,E15„4) , *N 

LF=5  (X,g (x) ) or  9 (Maxwellian); 

U is  used  if  LF=9;  U>0„ 

k)  NEPTH  (>2)  (15) 

No.  of  energy  points  for  0. 

l)  (E,0) , NEPTH  pairs,  (6E11„ 4) , *L*N 

Energies  increasing. 

m)  (Supply  only  if  LF=5;  skip  (m) 
and  (n)  if  LF=9) 

NPX  (>2)  (15) 

No.  of  X points  for  g(X) 

n)  (X,g (X) ) , NPX  pairs  (6E11.4) , *L*N 

Increasing  X. 

Repeat  h-1  (if  LF=9)  or  h-n  (if  LF=5) 
for  each  subsequent  distribution  within  a 
reaction  Repeat  items  e-n  for  each  sub- 
sequent reaction. 
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(X,I4,I5,I10,2E15*5),*N 


Item  4*8  - Type  8,  Cross  Section  for  Level  Excitation 

a)  KK,  NES,  ID,  EL,  EH 

NES=no„  of  energies  for  a,  , ;(>2)„ 

level  — 

b)  NLEVP  ( 1 <NLEVP  <^1 00 ) (15) 

No*  of  perturbed  levels: 

c)  (KLEV(K),  NEL(K),  (NLEVP  pairs)  (1615) 

KLEV=level  no*,  increasing;  not 

necessarily  consecutive  numbers* 

NEL=no«  of  energies  (from  highest 
down)  at  which  a level  is  given; 


2<NEL(K)<NES; 


d (NEL  (K)  ) 
dK 


<0* 


d)  E(1,2,„.,NES)  (7E11*4)  ,P*N 

Energies  for  cross  sections; 

decreasing  order*  Must  span  EH,  EL* 

e)  For  each  of  the  NLEVP  levels: 

0rf  (1,2,„„*,NEL(K))  (7E11 „4) , P*N 

Values  of  0 level  (barns)  at  E(l), 

E (2)  „ * .. E (NEL  (K)  ) . 

Start  a new  card  for  each  level*  If 
scaling  is  used,  first  card  for  each 
new  level  must  have  P*N  desired  for 
that  level* 

Item  4.9  - 'IVpe  0,  Angular  Distribution  of  Discrete-Level  Scattering 

a)  KK,  NES,  ID,  EL,  EH  (X, 14 , 15, 110, 2E15* 5) , »N 

NES=no.  of  energies  at  which  tables  of 

dp/d|j  are  supplied  (>2) 

b)  NLEVP  <1<NLEVP<1 00)  (15) 

No  of  levels  perturbed 
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c) 


(KLEV(K),  NEL(K),  (NLEVP  pairs) 


(1615) 


Item 


NEL(K)  is  no,  of  energies  at  which 
tables  of  (y,dp/dy)  are  supplied; 

KLEV's  increasing;  2<NEL  (K)<_NES 
d)  E(1,2,...,NES) 

Energies,  decreasing,  at  which 
distribution  tables  are  supplied, 
must  span  EH,  EL, 

For  each  perturbed  level,  supply  e.l,  e.2, 
repeated  for  NELlevel  energies: 

(e.l)  NMUP  (>2) 

No.  of  (y,dp/dy)  pairs 

(e.2)  (y , dp/dy) , (NMUP  pairs) 

First  and  last  y's  should  be  -1 
and  +1.  dp/dy  may  be  un-normalized. 

Repeat  package  (e.l,  e.2)  for  all  energies  E(l),.. 
for  this  level. 

Repeat  package  (e.l, e.2)  for  all  subsequent  levels 
4.10  - Type-10,  Source  Spectrum 

a)  KK,  NES 

Pert.,  No.  of  energies  (NES>2). 

b)  (E , S (E) ) , NES  pairs 

Energies  increasing.  Must  span 
problem  EELOW,  EEHIGH. 

S(l)  is  constant  source  density  from 
E (1)  to  E (2)  , etc. 

End  of  Input. 


(7E11.4) , P*N 


(15) 


(6E11.4) ,*L*N 


. , E (NEL (K)  ) 


(X,  14, 15) 


(6E11.4) ,*L*N 
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3.4  Tape  Utilization 


The  following  describes  the  function  of  each  tape  used  in  this  program 
Tape  numbers  refer  to  Fortran  logical  numbers.  All  tapes  are  used  in  the 
binary  mode. 

Tape  8 

The  processed  perturbation  data  tape  generated  by  program 
SAM IN  to  be  read  in  by  the  next  programs,  SAMSAM  and  SAMCAR, 
as  well  as  SAMGAM  for  a secondary  y problem. 

Tape  12 

A temporary  storage  tape  used  in  subroutine  FILE4P  for  con- 
struction of  CHI-table. 

Tape  14 

A temporary  storage  tape  used  in  subroutine  FILE5P  for  con- 
struction of  ENN-table. 

Tape  15 

A temporary  storage  tape  used  in  subroutine  FILE5P  for  con- 
struction of  ENN-table. 
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4.  PROGRAM  SAMSAM:  NEUTRON  TRANSPORT  PRE-PROCESSOR 


4 1 General  Discussion 

In  Chapter  2,  we  discussed  the  fact  that  SAMCEP  performs  Monte  Carlo 
samplings  from  sampling  densities  which  are  defined  to  be  least  upper  bounds 
(l.u  b ),  or  renormalized  l.u.b. , of  physical  sampling  densities  appropriate 
to  each  of  the  correlated  problems  considered 

Most  of  these  sampling  densities  are  generated  by  the  SAMCAR  code,  during 
the  course  of  Monte  Carlo,  only  as  the  need  arises. 

The  exceptions  are  angular  distribution  and  continuous  energy  distributions, 
for  which  the  sampling  tables  are  precomputed  in  the  relevant  energy  ranges 
prior  to  the  start  of  the  Monte  Carlo  calculation.. 

The  precalculation  of  these  sampling  tables  is  performed  by  the  SAMSAM 
code,  subsequent  to  the  generation  of  a BAND  tape. 

The  SAMSAM  code  accepts  as  input  the  unperturbed  data  in  the  form  of  a 
BAND  tape,  and  the  definition  of  perturbations,  as  specified  by  the  output  of 
the  SAMIN  program. 

The  sampling  tables  generated  by  SAMSAM  are  broken  up  into  energy  bands 
and  the  information  is  written  on  tape  12  in  a format  described  in  Appendix  C. 

A record  is  also  added  to  tape  8 (SAMINP  tape),  consisting  of  KEGEOM  and 
the  INBAND  array,  where  KEGEOM  is  the  length  of  the  largest  cross  section  data 
band  and  INBAND  is  of  dimension  NBAND . INBAND  (I)  = 1 or  0 if  any  tables  are 
present  in,  or  missing  from,  the  I—  band. 

4.2  Methodology 

Generally,  the  user  of  the  program  will  have  at  his  disposal  an  Element 
Data  Tape  (EDT,  output  tape  of  processor  SAM-X)  which  contains,  for  every  iso- 
tope in  the  problem,  a set  of  energy-dependent  interaction  cross  sections  The 
user  must  then  specify  each  of  the  material  compositions  appearing  in  the  prob- 
lem. A composition  is  defined  in  terms  of  atomic  concentrations,  , (in  units 
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24  3 

of  10  atoms/cm  ) of  each  isotope,  i,  in  the  composition.  These  may  be 

calculated  from  the  expression: 

-24 

= 10  x Avogadro's  number  x mass  density /atomic  weight.  For  com- 
pounds or  mixtures  the  concentration  of  each  component  must  be  specified. 

In  addition,  each  composition  must  be  identified  by  a composition  number. 

This  input  is  processed  in  conjunction  with  the  EDT  by  the  BAND  routine, 
which  generates  an  Organized  Data  Tape  (ODT) . The  ODT  contains  all  the  ele- 


mental cross  section  data  in  energy  band  structure  . During  the  tracking  pro- 
cess, the  data  retrieval  (DR)  routines  of  program  SAMCAR  will  use  this  informa- 
tion to  determine: 

1 The  probability  that  a particle  lias  an  interaction  in 
a region  of  given  composition, 

2 The  element  with  which  the  particle  interacts, 

3.  The  type  of  interaction  (absorption,  elastic  scattering, 
etc.)  occurring, 

4.  The  energy  and  direction  of  the  particle  after  interaction. 

The  sampling  CHI-  and  ENN-tables  of  a specific  element  are  tabulated  at 
the  energy  mesh  points  of  the  unperturbed  CHI-  and  ENN-tables  of  that  element. 

If  perturbation  type  2 (entire  element  perturbed)  is  present,  or  the  perturbation 
energy  ranges  overlap,  there  will  be  just  one  sampling  table  for  that  element 
in  the  affected  energy  band.  Only  if  there  is  an  energy  gap  in  the  basic 
nergy  mesh  with  no  perturbations  (type  2,  6 or  7) , will  there  be  separate 
;ampling  tables  for  the  separate  perturbation  energy  ranges  in  this  energy  band. 
At  each  energy  mesh  point,  a sampling  table  is  generated  by  subroutine  NVDOPE 
according  to  the  number  of  such  perturbations  affecting  this  energy.  Together 


The  energy  band  structure  may  bo  user  specified,  or,  optionally  computed  intoi  - 
nally  on  the  basis  of  available  memory  for  cross  section  storage  (i.e.,  auto- 
matic banding). 


4H 


with  the  sampling  tables  written  on  tape  are  the  perturbation  ID's  (i.e.. 


I 


IP's)  so  that  we  can  fill  in  the  location  of  sampling  tables  for  each  pertur- 
bation after  the  sampling  tables  are  read  into  memory  from  tape  during  the 
execution  of  program  SAMCAR. 

SAMSAM  first  calls  subroutine  BAND  to  process  the  unperturbed  element 
cross  section  data  into  energy  bands  and  to  store  the  processed  data  on  tape. 
Subsequently,  SAMSAM  reads  in  the  processed  data,  a band  at  a time,  into  the 
MASTER-array  (refer  to  organization  of  the  MASTER-array  in  Appendix  B)  and 
then  fills  the  perturbed  element  locations  into  MISTER-array  (refer  to  storage 
map  description  of  MISTER-array  in  Appendix  C) „ With  both  the  perturbed  and 
unperturbed  element  cross  section  data  in  memory,  SAMSAM  next  calls  subrout- 
ines SAMCHI  and  SAMENN  to  retrieve  all  perturbations  of  type  2,  6 or  7 in  the 
current  energy  band,  and  to  group  such  perturbations  of  the  same  element  to- 
gether For  each  element,  SAMCHI  (or  SAMENN)  will  call  subroutine  CHIA  (or 
ENNA)  to  interpolate  and  tabulate  the  unperturbed  and  perturbed  CHI-  (or  ENN-) 
distributions  in  the  energy  mesh  of  the  unperturbed  element's  CHI-  (or  ENN-) 
table.  CHIA  (or  ENNA)  in  turn  calls  subroutine  NVLOPE  to  generate  the  sampling 
CHI-  (or  ENN-)  table  at  the  given  energy  mesh.  The  length  of  the  sampling 
table  is  set  at  a constant  value  of  15  (including  0 and  1 in  the  case  of  CHI- 
table);  this  simplifies  the  logic  of  the  program,  and  the  chosen  length  is  con- 
sidered to  be  adequate  for  reasonable  superpositions  of  angular  or  energy  dis- 
tributions 

It  is  possible  that  perturbations  of  type  6 and  7 (angular  and  secondary 
energy  distribution)  exist,  yet  no  sampling  table  is  generated.  This  is  the 
case  when  the  perturbation  energy  range  (EH  and  EL)  falls  within  two  energy 
mesh  points  of  the  basic  energy  table. 
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4.3  Description  cf  Input  Data  for  Program  SAMSAM 
Item  1 - IODT  = 0 or  KEGEOM  (Format  15) 

IODT  = 0,  SAMSAM  will  call  subroutine  BAND  to  process 
element  data  tape  (output  tape  of  SAM-X)  ; 

IODT  = KEGEOM,  the  call  to  BAND  is  bypassed 
KEGEOM  is  the  total  length  for  all  the  element  data, 
as  computed  (and  displayed)  by  a call  to  BAND  in  a 
previous  execution  of  SAMSAM. 

★ 

Item  2 - NBAND  = total  number  of  energy  bands  supplied  (Format  15) 
or  NBAND  = 0 effects  the  automatic  banding  option 
The  user  supplied  value  determines  the  input  items  to  be 
supplied  next: 

If  NBAND  = 0,  supply  Items  3a  and  omit  3b 
If  NBAND >0 , omit  Items  3a  and  supply  3b 
Items  3a  - automatic  banding  input  (NBAND=0) 

Card  #1  - EXl , EXLAST  (Format  2E11.4) 

★ ★ 

high,  low  cross  section  energy  limits  ; 

Card  #2  - NROOM,  NEL  (Format  215) 

where  NROOM  is  the  number  of  words  of  memory  available 
for  the  longest  BAND  of  data  (internally  set  to  5000 
if  left  blank) ; and  NEL  is  the  number  of  different 
elements  in  the  present  calculation; 

Card  (s)  #3  - ID(I),  1=1, NEL  (Format  1615) 

complete  list  of  element  ID's  (ZZAAA)  in  any  order; 

* 

Currently  limted  to  19 

* * 

See  Appendix  H. 
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Item  3b  - user  specified  banding  (NBAND  >0) 

EBAND (I) ,1=1,  (NBAND+1);  energy  band  limits. 

(Format  7E11.4) 

Enter  the  energy  limits  (in  ev)  of  each  band  starting 

with  the  highest  energy  and  proceeding  to  the  lowest 

* 

energy  in  the  problem  „ Use  as  many  cards  as  necessary. 
End  of  data  if  IODT=KEGEOM. 

If  IODT=0,  continue  with  next  items. 

Item  4 - Composition  Identification  (Format  7110) 

NO  = Problem  number  - identification  of  present  run. 

NCOMP  = Composition  number  - total  number  of  compositions  in 
the  problem. 

NG  = 0 for  neutron  transport  problem. 

Item  5 - NE  - Number  of  elements  (Format  110) 

Enter  the  number  of  discrete  nuclides. 

Item  6 - Element  cards  (one  card  for  each  of  the  NE  elements) 
(Format  2110,  E15.6) 

IT  - blank 

ID  - an  integer  which  identifies  each  element. 

(Five  decimal  digits:  ZZAAA) 

Concentration  - enter  the  atomic  concentration 
24  3 

(10  atoms/cm  of  this  element  in  the  composition) 

Items  5 and  6 are  repeated  for  each  composition. 

End  of  data  for  program  SAMSAM. 


See  Appendix  H 


4 4 Tape  Utilization 

The  following  describes  the  function  of  each  tape  used  in  this  program 
Tape  numbers  refer  to  Fortran  logical  numbers.  All  tapes  are  used  in  the 
binary  mode 
Tape  8 

The  processed  perturbation  data  tape,  output  of  program  SAMIN. 

★ 

Data  on  this  tape  is  modified  and  expanded  by  SAMSAM  . This 
tape  will  be  read  in  by  next  program  SAMCAR,  and,  subsequently, 
by  SAMGAM,  if  a secondary  y problem  is  being  solved. 

Tape  9 

★ ★ 

A temporary  storage  tape  is  used  by  subroutine  BAND  if  NBAND  > 1. 

Tape  10 

The  organized  data  tape,  output  tape  of  BAND  with  cross  section 
data  in  band  structure. 

Tape  11 

The  neutron  element  data  tape,  output  tape  of  processor  SAM-X 
which  contains  a library  of  element  data  including,  at  least, 
all  elements  required  for  the  problem  being  run.  Subroutine 
BAND  reads  in  this  tape  to  process  its  data  into  energy  band 
structures . 

Tape  1 2 

This  tape  is  generated  by  program  SAMSAM.  It  contains  the 
sampling  CHI-  and/or  F.NN-table  in  the  same  band  structure  as 
the  unperturbed  data.  This  tape  will  be  read  in  by  the  next 
program  SAMCAR  if  perturbations  of  this  type  are  present. 

* 

Number  of  bands  and  energy  band  limits  are  appended  by  SAMSAM,  to  r>e  read  in 
by  SAMCAR  (maximum  number  bands  currently  19  ) 

★ ★ 

Although  NBAND  is  set  to  zero  (by  the  user)  to  invoke  automatic  bandinq, 
it  is  internally  reset  to  correspond  to  the  number  of  computed  BAND  intervals. 
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5.  PROGRAM  SAMCAR:  NEUTRON  AND  SECONDARY  GAMMA 


CORRELATED  MONTE  CARLO  TRANSPORT  CODE 
5.1  General  Discussion 

SAMCAR  is  the  main  correlated  Monte  Carlo  program  which  calculates  the 

★ 

transport  of  particles  through  matter  (both  time-dependent  or  time- indepen- 
dent) with  correlated  sampling  techniques.  It  is  composed  of  a series  of 
independent  routines  which  perform  the  following  three  basic  functions: 

1 Process  geometry  data  by  calling  GENI. 

2.  Process  the  unperturbed  and  perturbed  source 

** 

spectrum  by  calling  subroutine  SOUCAL 

3.  Perform  the  transport  calculations  and  score 
the  particle  fluxes  by  calling  subroutine  CARLO. 

SAMCAR  writes  the  particle  fluxes  by  supergroups 
of  all  problems  (unperturbed  and  perturbed)  by 
statistical  aggregates  on  tape  to  be  edited  by 
the  next  program  SAMOUT. 

Easically,  the  program  requires  as  input  a geometry  specification,  the 

elemental  composition  of  each  region,  and  a specification  of  the  location  and 

* * 

time-,  energy-,  angular  distributions  of  the  radiation  source  . The  program 
selects  individual  particles  from  the  given  source  distribution  and  tracks 
them  through  a series  of  interactions  within  the  geometry  until  such  time  as 
the  particle  history  is  terminated  The  tracking  of  a particle  car  be  terminated 
for  any  of  the  following  reasons: 


★ 

neutrons  and  secondary  gamma  rays 

* * 

for  primary  neutron  problems. 

* * * 

or  external  source  tape,  for  secondary  y problems. 
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1„  The  energy  of  the  particle  after  an  interaction 
falls  below  a specified  'cutoff  energy' 

2„  The  time  variable  of  the  particle  exceeds  a 
specified  “cutoff  time',, 

3„  The  particle  escapes  from  the  geometry  (crosses 
an  external  boundary) „ 

4„  The  particle  is  'killed' „ This  procedure  will  be 

explained  in  the  section  dealing  with  the  importance 

sampling  techniques  employed  in  the  program 

* 

For  each  scoring  region  traversed  by  a given  particle,  the  code  computes 
the  flux  per  unit  time  per  unit  energy  as  a function  of  energy  and  time  for 
each  problem^  The  flux  contribution  for  a given  particle  is  defined  as  its  ex- 
pected total  path  length  contribution  in  a region  divided  by  the  volume  of  the 
region,.  Individual  particle  flux  contributions  are  accumulated  so  that  the  end 
result  of  the  tracking  process  is  the  total  flux  in  each  region  in  a specified 
group  of  energy  and  time  bins„  At  the  user's  option  the  problem  can  be  made 
time  independent 

The  above  description  of  the  SAMCAR  program  is,  of  course,  a very  simpli- 
fied view  of  the  computational  procedure  The  following  sections  provide  a more 
detailed,  although  non-mathematical , description  of  each  part  of  the  computation.. 


* Refer  to  Sections  5„2„1  and  5„2  9 
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5 „ 2 Input  Requirements 


1 


5.2.1  Bounded  Point  Detectors;  Region  Terminology 

With  the  inclusion  of  bounded- f lux-at-a-point  (BFAP)  estimation,  the 

concept  of  a scoring  region  has  been  generalized  to  include  point  detectors. 

This  generalization  has  necessitated  an  expansion  of  "region"  terminology 

with  the  introduction  of  two  qualifiers:  physical  and  material.  Thus,  a 

point  detector  is  located  within  a "physical"  region;  its  flux  estimates  are 

* 

scored  in  an  input  designated  "scoring"  region;  and  its  sphere  of  influence 
is  determined  from  the  element  concentrations  in  its  designated  "material" 
region,  or,  as  a user  input  option,  given  explicitly 

For  each  detector,  the  user  specifies  a coordinate  location,  scoring 
region,  and  material  region.  The  detector's  physical  region  is  internally  de- 
termined from  its  position  coordinates.  If  the  material  region  is  unspecified, 
the  default  material  region  is  the  physical  region.  If  a negative  value  for 

it  it 

the  material  region  is  specified,  the  critical  radius  is  given  explicitly 
Note  that  any  physical  region  may  itself  be  designated  as  a scoring 
region,  and  the  element  concentrations  of  any  physical  region  (scoring,  non- 
scoring, or  even  the  escape  region,  which  may  comprise  a "phony"  material) 
can  serve  as  the  designated  material  region  for  any  detector 

In  the  sections  that  follow,  the  generic  term  "region"  will  refer  to  a 
physical  region,  where  the  specific  reference  is  clear  from  context.  Also, 
whenever  applicable,  references  to  scoring  regions  are  equally  valid  for  point 
detectors 


* 

Refer  to  Sec.  F.  1.1  of  Appendix  E. 

* * 
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Refer  to  Item  10  of  Sec.  5.3. 


5.2,2  Geometry  Input 


The  specification  of  the  geometry  and  its  input  format  is  described  in 
Appendix  A of  this  report,.  It  makes  use  of  the  Combinatorial  Geometry  Tech- 
nique „ The  geometry  data  must  be  prepared  in  units  of  centimeters  to  be  con- 
sistent with  the  cross  section  data , 

5,2 , 3 Importance  Sampling 
A General 

Importance  sampling  or  'weighting'  provides  the  user  with  a powerful 
method  of  controlling  the  direction  and/or  energy  of  particles  in  the  problem. 
The  purpose  of  a particular  problem,  for  example,  may  be  to  calculate  the 
fast-neutron  flux  in  a given  region  within  the  geometry.  Under  normal  circum- 
stances, the  probability  of  a source  neutron  reaching  that  region  at  high 
energy  may  be  quite  small,  requiring  a vast  number  of  source  neutrons  to  be 
tracked  before  an  adequate  statistical  estimate  of  the  flux  is  obtained. 

However,  with  proper  particle  weighting  the  code  can  be  made  to  concentrate  only 
n those  fast  neutrons  having  the  best  chance  of  reaching  the  chosen  region, 
conversely,  the  code  will  spend  little  time  tracking  neutrons  which  are  either 
traveling  in  the  wrong  direction  or  are  at  relatively  low  energy, 

Tii.  program  determines  the  relative  importance  of  a particle  from  a 
parameter  called  the  weight.  The  total  weight  (W)  of  a particle  is,  in  turn, 
determined  from  a combination  of  three  quantities  called  region  weight  (WR) , 

angular  weight  (W  ),  and  energy  weight  (W  ) , where  W = W x W x W . Values 
U t,  R fj  E 

of  w , W , and  W are  given  as  input.  The  following  brief  discussion  should 

R Si  E 

I i ov ilk.  the  user  with  a better  insight  into  how  these  weights  are  actually 

used  by  the  code. 

f.  i the  normal  uncorrelated  Monte  Carlo  problem,  a quantity  F is  assigned 
to  ■ ach  particle.  The  value  of  F is  1,0  for  a source  particle  The  code 
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calculates  the  probability  that  the  particle  will  reach  the  boundary  of  the 


r 


source  region  along  its  flight  path  without  collision.  This  value  is  called 
F".  The  probability  that  a collision  takes  place  in  the  region  is  then  1-F". 
However,  in  the  case  of  a correlated  problem,  we  would  have  instead  N problems 
(unperturbed  plus  (N-l)  perturbed  problems).  The  value  of  F for  each  problem 
is  normally  1.0  for  a source  particle.  But  when  the  source  spectrum  is  per- 
turbed (perturbation  type  10)  the  sampling  source  spectrum  is  the  envelope  of 
all  source  spectra  instead  of  any  given  source  spectrum.  F,  therefore,  is 
different  for  each  problem,  and  in  fact  F=1.0  for  the  problem  whose  source 
spectrum  happens  to  be  the  envelope  and  the  rest  will  have  a value  of  F<1  0 
which  is  the  ratio  of  a given  spectrum  over  the  envelope  spectrum. 

The  probability  that  a collision  takes  place  in  the  region  is  more  com- 
plicated in  the  case  of  a correlated  problem  when  the  macroscopic  total  cross 
section  is  perturbed  (perturbation  type  1,  2 and  3).  Here  we  do  not  sample 
from  any  particular  distribution  but  we  sample  from  the  envelope  of  all  distri- 
butions. The  probability  of  having  a collision  at  a distance  S for  the  ith 
problem  is  F^y^  e 'iSdS,  where  y^  = the  macroscopic  total  cross  section  of  the 

i problem,  and  the  probability  of  not  having  a collision  in  the  region  is 
-M.S1 

F"  = F^e  l , where  SI  is  the  distance  to  the  region  boundary.  The  sampling 

distribution  for  picking  a collision  position  is  the  envelope  of  all  F^y^e  ‘i 

* c 

for  i=l  to  N in  a given  region  . At  collision,  quantities  F^  are  calculated 
for  the  single  collision  event  for  each  of  the  problems  i,  as  the  ratio  of  the 
i<f'  collision  probability  distribution  function  F^y^e  yiS  over  the  single  samp- 

Q 

ling  distribution.  The  collision  mechanics  and  further  adjustments  to  the  F^’s 
will  be  discussed  elsewhere  (see  section  DR3) . In  any  case,  the  particle  coming 

Q 

out  of  collision  and  its  adjusted  F\’s  are  stored  for  later  processing  and 

* 

Bounded  estimation  procedure  modifies  this  usage  (refer  to  Sec  E.3.3 
of  Appendix  E. ) 
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do  not  affect  the  further  tracking  of  the  original  particle  The  further 
tracking  of  the  original  particle  consists  in  attempting  to  generate  further 
collisions  in  the  same  region,  and  when  they,  if  any,  are  exhausted,  to 
enter  into  the  next  region,. 

If  the  next  region  is  the  "escape  region",  the  tracking  of  the  particle 
terminates*  If  not,  the  entrance  into  the  new  region  will  modify  the  F^'s  by 
the  new  region  weight*  Suppose  that  the  source  particle  is  leaving  region  1 
where  the  weight  is  and  entering  region  2 where  the  weight  is  At  the 

boundary  the  ratio  of  weights  W^/W^  is  multiplied  by  FV's  and  the  particle  is 
given  a starting  value  of  F^=F"  W^/W^  in  region  2.  The  probability  of  the 
particle  reaching  the  next  boundary  of  region  2 uncollided  is  calculated  and 

multiplied  by  F^  for  each  problem  to  obtain  a new  value  of  FV„  Notice  that  if 

* 

W_,  is  large  compared  to  W^,  the  probability  is  high  that  no  latents  will  be 
produced  since  both  F^  and  F"  will  be  small  compared  to  unity*  In  fact,  a para- 
meter F,  is  an  input  to  the  program,  and  if  the  largest  of  these  F^'s,  on 
entry  into  a new  region,  is  lower  than  F , a random  number  between  zero  and  one 
is  picked  If  the  number  is  greater  than  the  largest  , the  history  is  ter- 
minated If  it  is  lower  than  the  largest  F^,  the  history  is  continued  with  all 
the  F.  s divided  by  the  largest  F.„  Thus,  by  establishing  weight  sets  properly, 
increased  numbers  of  collisions  can  be  forced  to  occur  in  important  regions, 
and  in  addition,  the  original  source  particles  will  continue  to  propagate 
through  the  geometry* 

Without  going  into  detail,  it  can  be  stated  that  a small  value  of  F 

z 

minimizes  the  number  of  kills  (increases  problem  running  time)  A large  value 
jximi  ei;  the  kills  (decreases  running  time  per  history)  but  increases  the 


Refer  to  Section  F in  5„ 2,3 


variance  (or  error)  of  the  answers,  requiring  more  source  particles  to  be  run 


The  optimum  value  of  F^  will  generally  lie  in  the  range  from  0-01  to  0 1 

In  order  to  facilitate  input  preparation,  the  three  components  of  the 
total  weight  will  now  be  discussed  separately,. 

Bp  Region  Weights 

A region  weight  (w  ) must  be  specified  for  every  region  in  the  problem 

Ordinarily,  these  weights  are  set  up  so  that  they  gradually  decrease  as  a 

particle  proceeds  from  the  source  toward  a region  in  which  the  flux  is  desired 

Weights  should  gradually  increase  in  regions  which  are  located  progressively 

further  from  the  'important'  regionso  On  the  input  forms  the  user  must  specify 

all  values  of  W to  be  used  in  the  problem,.  The  order  in  which  these  values 
R 

are  entered  determines  their  region  weight  number  (i*e,,,  the  first  value  of 

W is  labelled  weight  #1,  the  second  value  is  weight  #2,  etc  ) Then  for  each 
R 

region,  the  weight  number  to  be  used  in  that  region  must  be  specified 
C,  Angular  Weights 

By  using  angular  weighting,  it  is  possible  to  specify  preferred  directions 
for  a particle,  independent  of  the  region  location  of  the  particle.  The  user 
first  specifies  the  direction  cosines  (with  respect  to  the  coordinate  axes  of 
the  problem)  of  one  or  more  aiming  angles.,  These  vectors  serve  as  ' zero  direct- 
ions' about  which  angular  weights  will  be  given,.  Each  aiming  angle  is  assigned 
a number.  Next,  a set  of  angular  bins  is  specified  between  0°  and  180°,  with 
the  bin  boundaries  given  in  terms  of  their  cosines..  Thus,  if  one  desires  to 
specify  four  bins  of  equal  angle,  the  cosines  of  0°,  45°,  90°,  135°,  and  180° 
should  be  entered..  Then,  one  or  more  sets  of  angular  weight  values  are  given 
For  each  set,  a weight  value  (Wq)  is  specified  for  each  angular  bin,  Each  set 
is  also  assigned  a number,.  Finally,  for  each  region,  the  aiming  angle  number 
and  the  angular  weight  set  number  must  be  given*  To  illustrate  how  the  code 
uses  this  information,  assume  that  a given  region  has  been  assigned  aiming 
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angle  #1  and  angular  weight  set  #2.  A particle  enters  the  region  and  the 


ode  determines  that  the  particle  is  traveling  at  an  angle  ft  with  respect  to 
aiming  angle  #1  The  code  then  determines  which  angular  bin  encompasses  ft, 
goes  to  angular  weight  set  #2,  and  finds  the  value  of  W in  that  bin. 

Id  general,  as  the  particle  direction  (angular  bin)  becomes  more  impor- 
tant, the  value  of  w assigned  to  that  bin  should  decrease. 

D,  Energy  Weighting 

The  use  of  energ  ' weighting  enables  the  user  to  instruct  the  code  as  to 

which  particle  energies  are  most  important  in  a given  problem.  A set  of  energy 

tins  i l,  first  given,  where  the  bin  boundaries  are  listed  in  decreasing  order. 

then , one  or  more  energy  weight  sets  are  specified,  with  each  set  being  assigned 

i number  For  each  set  an  energy  weight  value  W must  be  given  for  each  energy 

E 

bin.  The  energy  weight  set  number  corresponding  to  each  region  is  then  given. 
Assume,  for  example,  that  a particle  of  energy  E is  in  a region  which  has  been 
assigned  weight  set  #1.  The  code  first  locates  the  energy  bin  which  encompasses 
E,  refers  to  weight  set  #1,  and  determines  the  value  of  W^  which  was  given  for 
that  bin  In  establishing  the  energy  weights,  the  more  important  energies 
i sld  generally  have  smaller  W values  than  the  less  important  energies 
Application  of  Weights  to  Tracking 

not  -d  earlier,  th&  total  particle  weight  is  the  product  of  W x W x W „ 

K 1(  b 

Me1  particle  weight  is  used  to  determine  the  number  of  collisions  that  a par- 
t i 1c  will  produce  given  that  it  has  a specified  energy  and  direction  in  a 
m-  phy  > al  region.  By  an  appropriate  choice  of  aiming  angle  and  angular 
•*  i jht  , part  icles  heading  downward  can  be  caused  to  have  more  collisions  than 
■ < 1*  heading  upward  in  the  same  region.  Thus,  more  computing  time  will 

be  pent  >n  the  ’important'  downward-directed  particles  and  their  descendents 
than  : Mie  ’less  important'  upward- directed  particles,, 


hO 


F„  Treatment  of  ’Latent'  Particles 


1 


If  a collision  does  occur,  the  program  calculates  the  energy  anc  dir- 
ection of  the  particle  emerging  from  the  collision,,  The  collided  particle 
is  stored  in  a latent  storage  table  together  with  the  N values  of  FC  (refer 

i 

to  section  A of  5„2„3),  and  will  be  picked  up  and  followed  as  though  it  were 
a source  particle  at  a later  time„  When  it  is  started  out  as  a real  particle, 
it  is  assigned  a new  F^  value  for  each  problem  which  is  the  product  of  F^  and 
the  ratio  of  the  weight  (product  of  region,  energy  and  angular  weights)  before 
and  after  collision,,  In  general,  the  new  F^'s  will  be  different  due  to  differ- 
ences in  energy  and  direction  of  flight  if  energy  and  angular  importances  are 
present,, 

The  program  stores  the  information  concerning  latents  in  a table  which  can 
hold  up  to  100  latents „ Prior  to  storage  a test  is  made  to  see  if  the  largest 
F^  of  the  latent  exceeds  the  input  value  of  If  so,  it  is  stored  If  not, 

a game  of  Russian  roulette  is  performed,  as  previously  discussed,  and  the  latent 
is  either  eliminated  or  has  its  FVs  renormalized  by  dividing  all  of  them  by 
the  largest  F^„ 

If  more  than  100  latents  are  generated  by  a source  particle,  the  program 
has  a 'squeeze'  routine  which  reduces  the  number  of  latents  in  a statistically 
valid  way. 

Although  the  use  of  importance  sampling  may  appear  to  be  a rather  compli- 
cated procedure,  the  user  will  generally  find  that  after  gaining  a little  ex- 
perience with  the  code  the  process  becomes  relatively  straightforward  and  easily 
applied  Certainly,  the  time  spent  in  learning  how  to  properly  apply  this  tech- 
nique will  be  well  worth  it  in  the  long  run,  since  it  enables  complex,  deep- 
penetration  problems  to  be  run  in  a reasonable  amount  of  machine  time 
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5.2,4  Source  Specification 


The  specification  of  the  initial  particle  source  provides  the  user  with 
several  options.  These  options  are  described  briefly  below, 

5,2.4  1 Internally  Generated  Primary  Neutron  Source 

Perturbation  of  the  source  is  allowed  only  in  changing  the  shape  of  the 
source  spectrum  itself,  and  the  input  of  the  unperturbed  and  perturbed  spectra 
must  be  in  the  form  of  histograms. 

A Spatial  Distribution 

Sources  may  be  generated  in  any  number  of  regions,  but  the  region  must 
be  single  bodies  and  restricted  to  SPH,  RCC,  BOX,  RPP,  TRC„  For  each  source  re- 
gion the  power  density'  (particles/volume)  must  be  given.  The  user  has  the 
option  of  normalizing  the  problem  to  a unit  source  or  to  the  total  input  power, 

B_ Angular  Distribution 

Sources  may  be  either  isotropic  or  monodirectional  but  the  same  angular 
distribution  must  be  used  in  all  source  regions,  (It  should  be  noted,  however, 
that  a source  may  be  generated  in  a finite  cone  by  specifying  an  isotropic  dis- 
tribution and  using  angular  weights  to  kill  particles  which  are  generated  outside 

the  desired  cone,) 

( ' Energy  Distribution 

The  source  energy  spectrum  is  specified  by  a histogram, 
source  energy  spectrum  contains  the  desired  energy  mesh  EVs 
the  histogram  in  interval  IP  to  E^  + ^ (i  e , a table  of  E vs, 
r Tinv  Distribution 

i it  :m  dependent  problem  is  to  be  run,  the  user  must 
m i i i the  integrated  source  up  to  each  time  (i  e , 

i ini  it  il  1 b*  deleted  from  time- independent  problems. 


f>2 

- - - ^ 


The  input  of  the 
and  the  height  of 
S(E)  is  required), 

supply  a table  of 
t vb/;  S(t)dt), 


5. 2.4.2  Secondary  Gamma  Source  from  Previously  Generated  Interaction  Tape 

Using  a neutron  interaction  tape  (see  Section  5.2.12),  and  gamma  ray 

production  data  supplied  by  SAM-X,  SAMGAM  can  generate,  internally,  sources 
of  secondary  gamma  radiation.  Note  that  the  interaction  tape  must  have  been 
generated  by  a precursor  primary  neutron  SAMCAR  calculation.  A detailed  de- 
scription of  Program  SAMGAM  is  given  in  Chapter  7. 

5.2  5 Time  Dependence 

SAMCEP  enables  the  user  to  compute  particle  fluxes  as  a function  of  time 
as  well  as  energy  and  position  The  user  selects  any  desired  time  bin  struc- 
ture for  the  problem  and  enters  the  bin  limits  in  consecutive  order  on  the  input 
forms,  starting  with  the  latest  time.  Output  fluxes  will  be  given  in  tnis  oin 
structure  in  the  edit.  During  the  tracking  process  the  code  computes  the  flight 
time  of  a particle  between  collision  points  from  its  velocity  (or  energy). 
Interactions  are  assumed  to  occur  instantaneously.  By  accumulating  the  flight 
times  for  each  particle,  the  code  is  capable  of  storing  particle  fluxes  in 
the  proper  output  time  bins. 

5.2.6  Output  Energy  Mesh;  Supergroups  and  Superbins 

During  tracking,  the  code  stores  fluxes  in  each  region  in  a set  of  energy 
output  bins  specified  by  the  user  The  number  and  width  of  these  bins  are  ar- 
bitraiy.  The  bin  limits  muo  be  given  consecutively  in  the  input,  starting  wi  > 
the  highest  energy.  The  upper  and  lower  bin  limits  must  be  preceded  by  minus 
signs.  The  reason  for  this  will  be  explained  shortly.  Care  should  be  taken 
to  insure  that  the  upper  energy  bin  limit  is  equal  to  or  greater  than  the  high- 
■ .t  source  energy  to  be  generated  in  the  problem.  A cutoff  energy  is  also  speci- 
fied, which  instructs  the  code  to  cease  tracking  any  particle  which  degrades 
!i  1 w this  energy.  The  user  should  be  certain  that  the  lowest  energy  bin  limit 
i equal  to  or  lower  than  the  cutoff  energy.  Tn  essence,  there  must  be  a bin 
availahli  to  store  every  possible  energy  in  ttie  problem. 


Some  calculations  may  require  more  computer  storage  than  is  available „ 

This  situation  can  be  alleviated  by  using  the  'superbin'  option  provided  by  the 
code,  This  option  divides  the  overall  output  energy  range  into  smaller  groups 
(called  supergroups)  and  the  cross  section  energy  range  into  bands „ The  inter- 
section of  band  and  supergroup  defines  a superbin„  The  code  then  treats  each 
of  these  superbins  separately,.  In  this  manner,  only  the  cross  section  data  for 
the  band  currently  being  treated  are  stored  in  the  memory;  only  fluxes  for  the 
supergroup  currently  in  memory  are  scored;  and  only  those  particles  having 
energies  in  the  corresponding  superbin  are  tracked,,  When  a particle  degrades 
to  a lower  bin,  its  parameters  are  stored  and  its  tracking  is  resumed  only  after 
all  higher  superbins  have  been  completed,,  The  output  supergroup  structure  is  de- 
fined by  the  user  by  placing  a minus  sign  before  those  output  energies  he  wishes 
to  designate  as  supergroup  limits,.  If  this  option  is  not  desired,  only  the  upper 
and  lower  energy  bin  limits  require  minus  signs,.  This  instructs  the  code  to 
treat  all  output  energies  as  part  of  a single  supergroup 
So  2, 7 Response  Functions 

SAMCEP  provides  response  function  options  which  allow  the  user  to  auto- 
matically transform  particle  fluxes  into  any  desired  flux-dependent  quantity 
(done,  heat  deposition,  etc„)  Assume,  for  example,  that  the  dose  is  required 
in  several  regions,.  The  user  supplies,  as  input  in  the  next  edit  program, 

SAMC'UT,  a flux-to-dose  conversion  factor  as  a function  of  energy.  For  each  re- 
gion, SAMOUT  multiplies  the  flux  ^ (E)  in  each  energy  bin  by  the  corresponding 
conversion  factor  D(E)  and  integrates  over  energy.  Thus, 
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mission  region  is  used  when  it  is  desired  to  run  a problem  in  two  steps.  This 
is  usually  done  for  very  deep  penetrations  or  for  unusual  geometric  configura- 
tions (such  as  ducts)  where  it  may  be  more  economical  to  run  the  problem  in 
stageso  The  designation  of  a transmission  region  is,  however,  optional , N>  ' 
that  the  program  is  capable  of  treating  up  to  10  different  transmission  reg 

An  escape  region  is  one  in  which  all  particles  that  enter  are  killed  It 
is  ordinarily  used  to  define  the  outer  limits  of  the  geometry  (i  e„,  the  com- 
plete geometry  is  enclosed  in  a large  region  which  is  designated  as  the  escape 
region) „ 

5,2.9  Scoring  Regions 

A scoring  region  is  one  in  which  the  flux  contribution  is  computed  for 

* * * 

each  particle  which  passes  through  it  , In  a non-scoring  region  no  such  com- 
putation is  made,  so  that  the  output  edit  provides  fluxes  only  in  those  regions 
designated  in  the  input  as  scoring  regions.. 

In  most  problems  it  is  desired  to  know  the  flux  in  every  region  separately, 
in  which  case  each  region  in  the  problem  would  be  defined  as  a scoring  region 
with  a different  number.  In  some  problems,  however,  two  or  more  regions  may  be 
completely  symmetric  with  respect  to  the  source,  in  which  case  the  fluxes  in 
those  symmetric  regions  could  be  combined  without  any  loss  of  information,  and 

* 

Not  applicable  for  point  detectors. 

* * 

See  Section  5C2012 

Recall  that  this  concept  has  been  generalized  to  include  point  detector*;. 

(See  Section  5.  2,  1) 


65 


in  fact,  an  improvement  in  the  accuracy  will  be  obtained,.  Each  of  these 


regions  then  would  be  designated  by  the  same  scoring  region  number,,  In  still 
other  problems  it  may  be  unnecessary  to  know  the  fluxes  in  certain  regions. 

These  should  then  be  given  scoring  region  number  zero,  which  tabs  them  as  non- 

scoring. 

Since  fluxes  are  only  stored  and  printed  out  for  scoring  regions,  it  is 
possible  to  reduce  both  the  size  of  the  edit  and  the  core  storage  requirements 
by  reducing  the  number  of  different  scoring  regions.  It  should  be  remembered, 
however,  that  once  a problem  is  run  it  is  impossible  to  recapture  any  flux  in- 
formation in  nonscoring  regions, 

5.2,10  Number  of  Histories  and  Statistical  Groups 

The  user  must  designate  the  total  number  of  source  particles  (histories) 
to  be  run  in  the  problem.  Although  a greater  number  of  histories  will  improve 
the  accuracy  of  the  answers,  it  will  also  increase  the  problem  running  time.  The 
user  must,  therefore,  strike  a balance  between  the  tolerable  errors  in  the 
answers  and  the  cost  of  running  the  problem.  In  complicated  problems  it  is 
usually  wise  to  run  a test  problem  of  100  to  200  histories  to  get  a 'feel'  for 
whether  particles  are  reaching  the  desired  regions.  If  they  are  not,  the  fault 
probably  lies  in  incorrect  importance  sampling  and  the  weights  should  be  ad- 
justed, If  the  test  problem  appears  to  have  run  'well',  then  the  number  of  his- 
tories can  be  increased  by  perhaps  a factor  of  about  10,  After  some  experience, 
the  user  can  generally  determine  the  correct  number  of  histories  to  run  in  a 
particular  problem 

In  running  the  problem  the  total  number  of  histories  is  divided  into  aggre- 
; ci t • . ..Hod  statistical  groups  This  is  done  in  order  to  compute  the  variance 
to)  standard  deviation)  of  the  fluxes.  All  particles  (and  their  latents)  within 
a group  are  tracked  before  another  qroup  of  particles  Fluxes  are  computed  and 
stored  on  tape  separately  for  each  group  The  size  of  the  statistical  group  is 
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constant  in  a given  problem  and  must  be  specified  in  the  input.  The  groui 
size  is  not  critical  but  should  be  small  compared  to  the  total  number  of 
histories  to  be  run.  About  20  statistical  groups  per  problem  generally  have 
been  found  to  be  adequate. 

5.2.11  Volume  Computations 

To  evaluate  the  flux,  the  track  length  in  a region  is  divided  by  the  vol- 

★ 

ume  of  the  region  . Provision  has  been  made  to  input  volumes  of  regions  if 
they  are  known.  It  often  happens,  however,  that  regions  described  by  the  Cor., 
binatorial  Geometry  technique  have  such  complex  shapes  that  an  analytic  volume 
computation  is  not  practical.  To  determine  the  volume  of  such  regions,  a 
routine  is  included  to  perform  a ray-tracing  numerical  integration  calculation 
of  the  volume. 

A point  X within  the  geometry  is  given  as  input.  Rays  are  fired  isotropi- 
cally from  the  point  X and  the  volume  of  each  region  is  then 


»•! 


<rH> 


where  and  R2  are  the  distances  along  the  ray  to  the  entrance  and  exit  con- 
tacts with  region  i.  'H'  is  the  total  number  of  rays  fired. 

5.2.12  Interaction  File 

SAMCAR  provides  the  user  with  a method  of  calculating  the  production  and 
transport  of  secondary  gamma  rays  arising  from  neutron  capture  or  inelastic 
scattering  events. 


The  flux  estimates  scored  at  point  detectors  are  dimensionally  equivalent 
to  track  length  per  unit  volume. 
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During  the  tracking  of  primary  source  neutrons,  all  interactions  which  are 


capable  of  producing  secondaries  can  be  stored,  as  an  option,  on  an  "interaction" 
file.  Among  the  stored  data  are  the  coordinates  of  the  collision  point,  the  energy 
and  weight  of  the  primary  neutron  multiplied  by  the  non-elastic  interaction  proba- 
bility, and  the  time  of  interaction.  Subsequently,  this  file  can  be  processed  in- 
ternally by  SAMGAM  and  be  converted  into  a source  of  secondary  gamma  radiation. 

The  interaction  file  can  contain  both  interaction  and  transmission  informa- 
tion The  14th  parameter  of  each  particle  record  indicates  the  type  of  event. 

5 . 3 Card  Input  Formats 

Definitions  of  all  input  quantities  and  corresponding  card  formats  are  given 
below  The  definitions  are  given  in  the  order  in  which  the  data  is  required  by 

the  code. 

Item  1 - Edit  Title  Card  (Format  20A4) 

This  title  will  be  printed  by  the  next  edit  program,  SAMOUT. 

1 1 en i 2 - Debug  Switches  SSDR1,  SSDR3  (Format  211) 

SSDRl  = 1,  debug  printout  from  subroutine  DR1 
= 0,  no  debug  printout 

SSDR3  = 1,  debug  printout  from  subroutine  DR3„ 

= 0,  no  debug  printout. 

I tern  j - Identification  Card  (Format  211,  3X,  A3,  18A4) 

I DBG  - 1,  debug  printout  for  tracking  of  particle. 

= 0,  no  debug  printout. 

* 

iODT  = 1 , do  not  call  BAND  . 

NAMK  identification  of  present  SAMCAR  run. 

k 

i'ii'  ;tput  of  a i ri'Viou  ly  run  BAND  in  SAMSAM  or  SAMGAM  must  be  available 
and  designated  as  Tape  10. 


Item  4 - LRN  = Random  Number  Initiator  (Format  5X,I15) 


If  blank,  the  default  option,  (=1) , will  be  used. 

Code  insures  that  LRN  is  positive,  odd,  and  satisfies 
the  condition:  (1 . LE. LRN.LE. 2 (48 ) -1) . 

Item  5 - LDUMMY  - an  arbitrary  number  (Format  IIP) 

Item  6 - GENI  (Geometry)  Input 

The  Combinatorial  Geometry  (CG)  input  discussed  in 
Appendix  A must  be  preceded  by  a header  card 
(Format  110,  E15.4,  A3,  13A4) . 

IPRINT  = 0,  print  out  body  and  region  data  which  follow 
= 1,  print  out  body  and  region  data  as  well  as 
the  internal  arrays  in  which  they  are  stored 


= 2,  suppress  all  geometry  printout 
SCALE  multiply  all  CG  dimensions  by  this  scale  factor 


NG 


Enter  0 for  a neutron  problem,  or  1 for 


a gamma  problem. 

Number  of  output  time  bins  (enter  0 for 
a time-independent  problem). 

Number  of  output  energy  bins. 

Number  of  flux  scoring  regions  plus  the 
number  of  detectors  (as  given  in  Item  9 below) . 

Number  of  distinct  region  weights. 

The  escape  region  number. 

The  first  transmission  region  number  (leave  blank 
l f no  transmission  regions  are  desired). 

Note  that  the  above  11  items  appear  on  a single  card,  the  first  three 
i:-'  Fomat  110  and  the  last  eight  are  Format  15. 

ransml  n-Intera  t ion  Information  (Format  1415) 


. n region  numbers, 
ne  or  if  no  trans- 
■ ! f ising  point  de- 


i . are  to  be  recorded  on 

o*.  a leave  blank 

1 ■ l,i  t !*  ■ ar  e to  be  recorded 

thi'iwi-.e  leave  blank 


NT 

NOUT 

NUMSC 

NRWL 

1REX 

IRT(l) 
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Item  9 - Number  of  Point  Detectors  (Format  110) 


NDET  - Number  of  point  detectors  for  which  bounded 
statistical  estimation  will  be  used.  If  no 
detectors  used,  leave  blank  and  omit  next  item. 
Item  10  - Detector  Coordinates,  Scoring  Region,  "Material11 
Region,  and  (optionally)  critical  radius. 

(Format  3E14.6,  215,  E14.6) 

Use  one  card  for  each  of  the  NDET  detectors  (NDET>0). 


XAD  ' 

YAD 

ZAD 


The  X,  Y,  Z coordinates  of  the  detector 


ISCA  - The  scoring  region  of  the  detector.  A blank 
entry  will  suppress  scoring  for  the  detector. 
Otherwise,  any  number  in  the  range  (1,NUMSC), 
where  NUMSC  is  specified  in  Item  7,  is  valid, 
but  normally  given  as  (NUMSC-NDET+1)  through 
(NUMSC) . 

IRDET-  The  detector  "material"  region  (if  given  >0) . 

The  radius  of  the  detector's  sphere  of  influence, 
i.e  , "critical  radius",  is  computed  internally 
as  the  reciprocal  of  the  sum  of  the  concentrations 
of  the  elements  in  this  specified  region  A blank 
entry  will  assign  the  region  in  which  the  detector 
is  located,  as  specified  by  its  coordinates  (i.e  , 
its  physical  region).  A negative  entry  implies  cri- 
tical radius  specified  by  user  as  next  entry  on  card 
CRAD  - Critical  radius  (ignored  if  IRDET'-O)  . 
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Item  11  - Cutoff  Information  (Format  2E14.5,  14X,  2E14*5) 


ECUT  Low  energy  cutoff  (ev)»  Tracking  of  a particle 

* 

is  terminated  if  its  energy  degrades  below  ECUT  „ 

ETH  Thermal  energy  if  a thermal  group  is  required 

ETH  must  be  within  the  energy  limits  of  the  problem* 

FZ  See  discussion  in  Section  5.2.3. 

EHIGH  High  energy  cutoff  (ev) . This  should  be  less  than 
or  equal  to  the  highest  energy  for  which  cross  sec- 

it 

tions  are  available  o 

Item  12  - Output  Energy  Bins  (Format  5E14„5) 

These  cards  give  the  boundaries  of  the  output  energy  bins  (ev)  used 
* 

for  the  flux  edit  . There  should  be  five  entries  per  card  with  a total 
of  (NOUT+l)  entries  The  energies  should  be  listed  from  high  to  low 
with  the  lowest  energy  equal  to  or  less  than  ECUT.  The  first  and  last 
entries  should  be  negative*  If  supergroups  are  used,  any  number  of  inter- 
mediate energies  may  also  be  negative*  The  absolute  value  of  the  high- 
est output  bin  boundary  must  be  >i:HIGH. 
item  13  - Output  Time  Bins  (Format  5E14.5) 

These  cards  give  the  boundaries  of  the  output  time  bins*  There  should 
be  five  entries  per  card  with  a total  of  (NT+1)  entries*  However,  if 
NT=0,  omit  Item  13  entirely*  Times  should  be  entered  from  high  to  low 
and  the  last  entry  must  equal  0*  The  first  entry  defines  TCUT,  the  time 

cutoff* 

1 1 ■ ^ 1 4 - Region  Weights  (Format  5E14*5) 

These  cards  give  all  of  the  region  weights  needed  in  the  problem*  They 
art  entered  five  to  a card  with  a total  of  NRWL  entries.  The  weights 
need  not  be  entered  monotonically  by  value,  but  their  order  determines 
the  region  weight  numbers  (i  e„  entry  one  is  weight  #1,  etc*). 


Si  i Appendix  II 
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Item  15  - Region  Specifications  (Format  615) 

Use  one  card  per  region  with  a total  of  NRMAX  cards.  The  first 
card  applies  to  region  1,  the  second  to  region  2,  etc. 

ISC  Scoring  region  number  in  which  the  fluxes 

in  this  physical  region  are  to  be  stored. 

Several  regions  may  be  assigned  the  same  ISC 
number.  If  ISC=0,  fluxes  will  not  be  scored. 

NREG  Number  of  the  composition  to  be  found  in  this 

region. 

IRW  Region  weight  number  assigned  to  this  region. 

A weight  number  is  given  by  its  position  in  the 
list  of  region  weights. 

IEW  Energy  weight  set  number  assigned  to  this  region. 

If  IEW=0,  there  is  no  energy  weighting  in  this 


region. 

IAN  Aiming  angle  number  assigned  to  this  region  If 

IAN=0,  there  is  no  angular  weighting  in  this  region. 

IANG  Angular  weight  set  number  assigned  to  this  region. 

Item  16  - Energy  Weight  Specification  (Format  2110) 

NEWL  Number  of  energy  bins  for  energy  weighting. 

NEW  Number  of  distinct  energy  weight  sets.  If  NEWL  and 

NEW=0,  the  problem  contains  no  energy  weighting  and 

Items  17  and  18  are  omitted. 

Item  17  - Bin  Limits  for  Energy  Weights  (Format  5E14.5) 

Enter  the  boundaries  (ev)  of  the  energy  bins  to  be  used  for  energy 
* 

weighting  . There  should  be  five  entries  per  card  with  a total  of 
(NEWL+1)  entries.  The  energies  should  be  entered  in  decreasing  order. 
The  lowest  bin  limit  should  be  less  than  or  equal  to  EOUT. 


Appendix  H 
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Item  18  - Energy  Weight  Sets  (Format  5E14.5) 


The  energy  weight  value  in  each  of  the  above  energy  bins  should 
be  entered.  One  or  more  sets  of  energy  weights  may  be  entered. 

Each  set  should  contain  NEWL  entries  and  a new  card  should  be  used 
to  start  each  set.  There  should  be  a total  of  NEW  sets.  The  order 
in  which  the  sets  are  entered  determines  the  energy  weight  set 
numbers  (the  first  set  is  weight  set  #1,  etc.). 

NOTE  - Omit  items  17  and  18  if  no  energy  weighting. 

Item  19  - Angular  Weight  Specif ications  (Format  3110) 

NAIML  Number  of  distinct  aiming  angles. 

NUMANL  Number  of  angular  bins  for  angular  weighting. 

NUMANG  Number  of  distinct  angular  weight  sets. 

If  the  problem  contains  no  angular  weighting,  enter  0 for  the  above 
three  quantities  and  omit  Items  20,  21  and  22„ 

Item  20  - Aiming  Directions  (Format  3E14.5) 

Enter  the  direction  cosines  of  each  aiming  angle  with  respect  to  the 
X,  Y,  Z coordinates.  Use  a total  of  NAIML  cards. 

Item  21  - Bin  Limits  for  Angular  Weights  (Format  5E14.5) 

Enter  the  boundaries  of  the  angular  bins  to  be  used  for  angular 
weighting.  Boundaries  are  given  in  terms  of  the  cosines  of  the  angles 
between  aiming  direction  and  particle  direction,  with  the  first  entry 
equal  to  1.0  and  the  last  entry  equal  to  -1.0.  There  would  be  a total 
.f  (NUMANL+1)  entries. 
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Item  22  - Angular  Weight  Sets  (Format  5E14.5) 


The  angular  weight  value  in  each  of  the  above  angular  bins  should 
be  entered.  One  or  more  sets  of  angular  weights  should  be  entered. 
Each  set  should  contain  NUMANL  entries  and  a new  card  should  be 
used  to  start  each  set.  There  should  be  a total  of  NUMANG  sets 
The  order  in  which  the  sets  are  entered  determines  the  angular 
weight  set  numbers. 

NOTE  - Omit  Items  20,  21  and  22  if  no  angular  weighting. 

Item  23  - Source  Specifications  (Format  3110) 

NSR  Number  of  different  source  regions  in  the 

problem.  If  NSR=0,  an  external  source  tape 
is  used  and  no  further  source  input  is  re- 
quired (i.e.,  omit  Items  24-28). 

IFLAG  Number  of  energies  used  to  define  the  source 
spectrum. 

ISW  If  ISW=0 , fluxes  will  normalize  to  one  source 

particle.  If  ISW=1,  fluxes  will  be  lormalized 
to  the  total  source  power  as  implie  by  Item  24. 

Item  24  - Source  Regions  (Format  IIP,  E20.8,  IIP) 

One  card  is  required  for  each  source  region  with  a total  of  NSR 
such  cards. 


ISR  Geometrical  region  number. 

P Power  density  in  the  region  (source 

particles  per  unit  volume) 

ISO  If  ISO=0,  the  source  will  be  emitted  isotropically. 

If  ISO=l,  the  source  will  be  monodirectional , with 
the  direction  specified  on  Item  28.  The  value  of 
ISO  i ' or  1)  must  be  the  same  for  all  source  reqions. 
NOTE:  Only  single  body  source  reqions  are  allowed.  Moreover, 

the  choice  of  bodies  is  restricted  to  the  following: 

SPH,  RCC,  BOX,  RPP. 
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Item  29  - Volume  Computation  Parameters  (Format  3110,  3E10  2) 


NMA  The  number  of  rays  to  be  fir  ad  for  volume 

computation. 

NST  The  number  of  rays  in  a statistical  aggregate 

for  the  purpose  of  error  computation 
IRSAVE  The  region  number  from  which  rays  are  to  be  fired. 

X, Y , Z The  X, Y, Z coordinates  of  the  point  from  which 

rays  are  to  be  fired 

NOTE  - if  volumes  are  precomputed  and  no  volume  computation  is 
needed.  Item  29  is  a blank  card. 

Item  30  - Precomputed  Volumes  (Format  I10,E10.B) 

IP  Region  number  of  precomputed  volume. 

VP  Precomputed  Volume. 

NOTE  - if  no  precomputed  volumes  are  supplied  Item  30  is  a blank 
card  As  many  precomputed  volumes  as  desired  (in  ascending  region 
number)  may  be  supplied.  In  any  case,  the  last  card  must  be  blank. 

If  no  volume  computation  and  no  precomputed  volumes  are  desired, 
then  2 blank  cards  corresponding  to  Items  29  and  30  are  supplied. 

This  effects  the  default  option  of  1.0  for  all  the  volumes. 

Item  31  - Fast  Slowing  Down  Option 

A blank  card  will  cause  this  option  to  be  ignored.  Otherwise,  for  air 
over  ground  problems  only: 

(a)  (Format  4E10.4) 

ECTE  Neutron  energy  (eV)  below  which  option  is  invoked. 

STMAX  Lateral  extent  in  cm.  (from  source  center)  of  slowing 
down  region  (internally  squared). 

ZS  Source  height  in  cm. 

RS  Source  radius  in  cm.  (assumed  spherical  source  region) . 


J 
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(b)  (Format  I5/(6E10.3)) 


NALT  Number  of  air  layers  of  constant  density  plus  one 

(ground  layer) 

[ZBD(K) ,DFR(K) ,K=1,NALT]  - 

Pairs  of  values  for  lower  layer  boundary  and  corres- 
ponding layer  density  factor  (ground  and  contiguous 
air  layer  factors  set  to  lo0  internally). 

ZBD(NALTtl)  Upper  layer  boundary  of  uppermost  air  layer,  defining 
upper  boundary  of  slowing  down  region. 

End  of  data. 

5.4  Organization  of  the  SAMCAR  Program 

The  SAMCAR  program  comprises  a small  main  program,  a supervisory  routine 
(the  effective  main  routine) , and  associated  subroutines  The  main  program 
(SAMCAR)  establishes  certain  array  dimensions  and  initializes  several  flags. 
The  supervisory  routine  (AMCAR)  controls  the  entire  calculation,  calling  the 
input  processing  routines  (GENI  and  INPUTD) , and  the  principal  Monte  Carlo 
routine  (CARLO).  In  addition,  the  supervisory  routine  controls  the  operation 
of  superbins  A number  of  labeled  common  blocks  are  used,  some  of  which  are 
described  in  detail  in  subsequent  sections  of  this  report. 
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5.4.1  Description  of  Routines 

This  section  gives  a brief  description  of  every  routine  in  the  program. 
The  routines  which  have  an  asterisk  (*)  superscript  are  of  major  importance. 
Routines  which  follow  the  main  program  description  are  given  in  alphabetical 
order.  The  brief  descriptions  are  followed  by  detailed  descriptions  of  the 
important  routines. 


ROUTINE 

SAMCAR 


AMCAR* 


ANISOT 


ARG* 


DESCRIPTION 

The  main  program.  It  establishes  the  dimensions  of  the 
MASTER  and  MISTER  arrays,  initializes  certain  flags,  and 
reads  the  first  binary  record  of  information  from  the  pro- 
cessed perturbation  data  tape  (created  by  SAMIN  and  modi- 
fied by  SAMSAM  and  SAMGAM  for  a secondary  gamma  problem) . 
Control  is  then  passed  to  the  supervisory  routine,  AMCAR. 
The  control  routine  for  the  Monte  Carlo  calculations. 

This  routine  reads  the  remainder  of  information  from  tapes 
written  by  the  previous  jobs  SAMIN  and  SAMSAM  or  SAMGAM 
for  a secondary  gamma  problem.  It  manipulates  the  rgan- 
ized  data  tape  (ODT)  and  the  supergroup  tapes.  This  rou- 
tine also  writes  on  tape  the  scores  for  each  statistical 
aggregate  for  later  use  by  the  edit  program  SAMOUT 
A routine  to  pick  a CHI  value  which  will  give  the  cosine  of 
scattering  angle,  0,  in  the  case  of  anisotropic  scattering. 
The  general  angle  reselection  routine  of  the  bounded 
estimation  procedure. 
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ROUTINE 


DESCRIPTION 


ARP REP 


ASSIGN 


CARLO* 


CARSCA* 

DIREC 


A subordinate  angle  reselection  routine  called  by  ARG. 

This  routine  implements  the  special  reselection  procedure 
which  is  required  for  "hydrogen- type"  neutron  scattering, 
i.e.,  scattering  in  which  the  effective  mass  of  the  re- 
coiling nucleus  is  <1.0.  The  relevant  analysis  is  given 
in  Appendix  E. 

A routine  to  assign  IP's  (perturbations)  of  specific 
types  of  corresponding  problems.  It  returns  an  array 
IPROB (I , J) , where  1=1,  total  number  of  problems;  and  J=l, 

11.  IPROB (I,11)=NIP,  number  of  IP's  affecting  the  Ith 
problem.  IPROB (1,1  to  NIP)  gives  all  the  IP's  for  the 
Ith  problem. 

The  correlated  Monte  Carlo  calculation  routine.  The 
routine  performs  tracking,  and  controls  importance  sampling, 
collision  events  (via  calls  to  CARSCA) , and  track  length 
scoring. 

This  routine,  called  by  CARLO,  controls  collision  events. 
Using  a direction  vector  W as  a polar  axis,  a dire-tion  W' 
is  generated,  the  cosine  of  whose  polar  angle  is  CSTHT  and 
with  random  azimuthal  angle.  Thus,  new  direction  cosines 
W'  are  computed  such  that 
W*W'=CSTHT 

The  routine  is  used  in  computing  the  new  direction  after 
a scattering. 


BO 


ROUTINE 


DESCRIPTION 


DR1*  This  routine  calculates  the  unperturbed  and  perturbed 

macroscopic  total  cross  sections  as  a function  of  energy 
for  a given  composition.  It  also  computes  the  sampling 
distribution  for  picking  the  element  of  reaction. 

DR3*  The  collision  mechanics  routine.  This  routine  picks  the 

element  of  reaction  and  the  reaction  type.  For  Klein- 
Nishina  scattering,  it  governs  flux  estimation  at  point 
detectors  via  calls  to  FLUP. 

DR31*  This  routine  is  called  by  DR3  when  an  elastic  scattering 

event  has  been  picked.  It  governs  flux  estimation  at  point 
detectors  via  calls  to  FLUP,  determines  energy  and  angle 
after  scattering,  and  governs  the  weight  adjustment  of  each 
problem  (unperturbed  and  perturbed) . 

DR32*  This  routine  is  called  by  DR3  when  an  inelastic  event  has 

been  picked.  It  performs  the  same  functions  as  DR31. 

DR33  This  routine  computes  the  probability  of  scattering  ly  a 

given  angle  in  the  lab  system  and  the  resulting  energy.  It 
calculates  the  common  factor  (lab/c. m.  Jacobian)  and  all 
DR35  or  DR36  for  the  problem  dependent  c.m.  factor,  fo 

* 

elastic  or  discrete  inelastic  scattering,  respectively  . 

It  is  utilized  by  FLUP  in  the  BFAP  estimation  procedure, 
and  by  DR34  in  the  angle  reselection  procedure. 


Continuum  inelastic  is  currently  treated  as  isotropic  in  the  c.m. 
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ROUTINE 


DESCRIPTION 


DR34 


DR35 


DR36 


FLUP* 


GARB 


GENI 


GETIR 


This  routine,  called  by  DR31  and  DR32,  initiates  the 
angular  reselection  procedure  by  performing  the  prelimin- 
ary position  checking  (see  Appendix  F)  and  then  calling  ARG 
and  DR33  for  reselection  and  weight  adjustment. 

This  routine,  called  by  DR33,  SAMPIC  and  XWADJ,  computes 
the  weight  adjustment  factor  for  every  problem  (PAIPP  array) 
given  a value  of  CHI  for  an  elastic  scattering  event. 

This  routine,  called  by  DR33  and  XWADJ,  is  the  inelastic 
scattering  counterpart  of  DR35. 

This  is  the  f lux-at-a-point  estimation  routine.  It  is  called 
by  AMCAR  for  source  particles  and  "FLUP"  latents,  and  by 
DR3,  DR31,  and  DR32  for  Klein-Nishina , elastic,  and  in- 
elastic collisions,  respectively. 

This  is  one  of  the  routines  of  the  Combinatorial  Geometry 
Package  (CGP).  It  is  subordinate  to  the  principal  CGP  input 
processor,  GENI,  and  is  utilized  whenever  the  ARB  body  is 
specified  on  input. 

The  major  geometry  input  processing  routine.  The  routine 
reads  geometry  data,  checks  for  errors,  and  puts  the  data 
into  the  MA  and  FPD  arrays  in  the  form  required  by  the 
tracking  routines. 

Another  member  of  the  CGP  (as  signaled  by  its  leading 
character).  This  routine,  called  by  INPUTD,  determines 
the  physical  region  of  a point  detector,  given  its  position 

coordinates. 
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' 45$  .#  Vi' 


ROUTINE 


DESCRIPTION 


GG 


GP 


G1 


INGRAL 

INPUTD 


PACK* 


PICK* 


The  distance  calculating  routine  of  the  CGP.  Given  a 
position  X,  a direction  W and  a body  number,  the  routine 
computes  the  two  distances  RIN,  ROUT  measured  from  X to 
the  body. 

An  auxiliary  routine  called,  as  an  option,  by  G1  for 
debugging  printout. 

The  main  geometry  tracking  routine.  Given  a position  X, 
a direction  W,  and  a region  IR,  the  routine  will  calculate 
the  distance  'S'  from  the  point  X to  the  next  region  in 
the  direction  W.  The  routine  also  determines  IR',  the 
next  region  to  be  encountered. 

An  auxiliary  routine  called  by  SOUPIC  to  normalize  histograir 
The  Monte  Carlo  input  routine.  The  routine  reads  Monte 
Carlo  input  data  and  stores  it  in  the  MASTER  array  for  use 
by  the  calculation  routines.  Region  and  point  detector 
specification,  importance  sampling  data,  and  output  energy 
meshes  are  handled  by  this  routine. 

A routine  to  store  the  position,  direction,  energy  and 
weight  data  for  a particle  into  24  computer  words.  The 
24-word  records  are  used  for  latent  storage  and  for  output 
onto  the  interaction/transmission  tape. 

The  control  routine  for  supergroup  latents.  This  routint 
stores  particles  whose  energies  are  not  included  in  the 
supergroup  currently  being  processed.  These  particles  are 
stored  by  PICK  in  either  the  central  memory  or  on  tapes  de- 
pending on  the  amount  of  core  available  for  the  particular 
problem. 


H i 


ROUTINE 


DESCRIPTION 


RANF 

MPIC 


SAMP IE 


SEEK* 

jUCAL* 


soupic* 


The  random  number  generator. 

A routine  to  pick  a CHI  value  (which  will  give  the  angle  of 
scattering  at  a given  energy)  from  the  sampling  CHI-table. 
This  routine  (via  DR35)  adjusts  the  perturbation  weight  of 
each  problem  by  multiplying  into  each  problem  weight  the 
ratio  of  probability  density  of  selecting  this  CHI-value 
from  each  particular  distribution  over  the  sampling  distri- 
bution. 

A routine  to  pick  an  ENN-value  (the  emerging  energy  in  the 
event  of  continuum  inelastic  scattering)  from  the  sampling 
ENN-table.  This  routine  also  adjusts  the  perturbation  weight 
of  each  problem  by  multiplying  into  each  problem  weight  the 
ratio  of  probability  density  of  selecting  this  ENN-value  for 
each  particular  distribution  over  the  sampling  distribution. 
Given  a vector  ET  with  elements  monotonically  decreasing  or 
increasing,  and  a variable  E,  this  routine  will  search 
through  the  ET  and  determine  the  bin  containing  E. 

The  input  processor  for  the  source  information  required  by 
the  Monte  Carlo  routines  for  primary  neutron  problems.  The 
routine  processes  energy  and  time  spectrum  data,  and  source 
region  data.  The  data  are  stored  in  the  MASTER  array. 

A routine  to  generate  initial  source  particles.  The  routine 
uses  the  source  information  processed  by  SOUCAL  and  delivers 
position,  energy,  time,  direction,  and  weight  data  f >r  source 
particles.  It  also  roads  in  source  data  from  an  external 
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source  tape. 


ROUTINE 


DESCRIPTION 


TALLY* 


TERP 

TIME 

TRALA* 

TROPIC 

UNPACK* 

UNPR 

VCALC 

WR'1’14 

XWADJ (I EL) 


A summary  routine.  The  routine  prints  a one  line  summar 
of  results  for  each  statistical  aggregate.  Quantities  s 
as  number  of  collisions,  absorptions,  degradations,  birt 
and  deaths  are  printed  for  each  aggregate. 

A linear  interpolation  routine. 

This  routine  interrogates  the  machine  clock  routine.  It 
is  utilized  by  AMCAR,  TALLY,  and  VCALC. 

This  routine  is  called  by  FLUP  to  compute  the  problem  de- 
pendent array  of  optical  pathlengths  between  collision  ( 
source)  point  and  a detector  point. 

A routine  to  generate  a vector  of  direction  cosines  from 
isotropic  distribution. 

A routine  to  unpack  a 24-word  vector  containing  position 
direction,  energy,  .d  weight  information. 

A routine  to  retrieve  from  the  MASTER  array  six  integer 
variables  and  to  deliver  the  six  integer  variables  store 
in  the  common  block  labeled  REGPAR. 

The  volume  computation  routine.  Region  volumes  are  com- 
puted by  numerical  integration. 

A routine  to  write  24-word  records  onto  tape.  The  routi 
is  called  whenever  a transmission  or  interaction  is  to  1 
put  on  tape. 

This  routine  (in  conjunction  with  DR35  and  DR36)  adjusts 
perturbation  weight  of  each  problem  in  case  a perturbat! 
of  angular  distribution  exists.  IEL=1,  for  elastic  sca< 
ing;  IKL=2,  for  inelastic  discrete  scattering.  This  roi 
tine  is  called  only  when  there  is  perturbation  but  no 
samp] ing  table  is  available.  The  selection  of  angle  of 
scattering  in  this  case  is  from  the  unperturbed  distrib 
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; ^tive  Main  Program  - Subroutine  AMCAR 


AMCAR  controls  the  Monte  Carlo  subroutine  (CARLO)  and  arranges  storage 
'location  for  supergroups.  As  noted  in  Section  4.2,  BAND  arranges  cross 
sections  in  certain  energy  bands.  The  output  energy  bins  are  also  arranged 
in  certain  output  supergroups.  Cross  section  input  corresponding  to  a sing 
band  can  be  stored  in  computer  memory  at  a given  time.  Scores  correspondin' 
to  a single  output  supergroup  can  be  stored  in  computer  memory  at  any  given 
time.  The  two  meshes  (bands  and  output  supergroups)  do  not  necessarily  co- 
incide. A combined  mesh  defines  a set  of  superbins. 

AMCAR  starts  by  reading  in  the  highest  energy  cross  section  band,  and  ] 
arranging  the  memory  for  the  highest  energy  output  supergroup.  The  highest 
of  the  low-energy  bounds  of  these  energy  ranges  defines  EBL,  the  low  energy 
the  superbin  currently  treated. 

The  program  then  calls  the  source-picking  routine  SOUPIC,  and  examines 

the  energy  E.  If  E_<EBL,  the  particle  is  stored  as  a latent  by  calling  the 

t Mi-j  part  of  subroutine  PICK,  and  SOUPIC  is  called  again.  When  E> EBL,  s 

* 

i utine  CARLO  is  called  . CARLO  tracks  the  particle,  scores  contributions 
f luxes  when  needed,  and  distributes  collisions  along  their  track.  Particle 
c'i  i out  of  collision  are  also  tracked  if  their  energy  is  above  EBL;  part 
ig  out  of  collision  with  E<EBL  are  stored  as  latents  by  calling  PICK. 

1 is  finally  returned  to  AMCAR,  which  proceeds  to  the  next  source  partis 
This  process  is  repeated  until  a complete  statistical  aggregate  of  particle 
i boon  treated  for  the  highest  energy  superbin. 


If  the  run  involves  point  detectors  (i.e.,  NDET/O) , the  call  to  CARLi 1 
is  preceded  by  a call  to  FLUP. 


At  this  point,  AMCAR  switches  to  the  next  energy  range  by  either  reading 
a new  band  of  cross  section  data,  or  by  writing  out  on  tape  the  set  of  scores 
obtained  and  preparing  the  memory  layout  for  the  next  supergroup,  or  both. 

It  then  proceeds  to  call  the  retrieval  section  of  subroutine  PICK,  which  picks 

latents  from  previous  superbins.  If  E<_EBL , the  particle  is  stored  again  as  a 

* 

latent  by  calling  PICK.  If  E>EBL,  CARLO  is  called  . The  procedure  continues 
until  all  latents  have  been  examined,  at  which  point  AMCAR  switches  to  the 
next  superbin,  etc.,  until  the  low-energy  cutoff  is  encountered.  When  this 
occurs,  AMCAR  switches  to  the  highest  superbin,  and  proceeds  to  treat  the  next 
statistical  aggregate  of  particles.  The  calculation  terminates  when  a history 
number  exceeds  the  cutoff  value  NHIST  specified  on  input.  A 'blank'  interaction 
record,  with  NHIST=NHIST+1 , is  written  on  the  interaction  tape  and  all  tapes 
are  rewound. 

Subroutine  ARG 

This  is  the  general  angle  reselection  routine  of  the  bounded  flux-at-a- 
point  (BFAP)  estimation  procedure.  Together  with  its  auxiliary  subroutine 
ARPREP,  it  implements  the  algorithms  described  in  Appendices  E and  F.  A non- 
mathematical  description  is  given  in  this  section. 

Subroutine  ARG  is  called  by  SOUPIC  when  a source  direction  has  been  pi  :ked, 
or  by  DR34,  following  the  selection  of  a post-scattering  direction.  In  either 
case,  this  initially  selected  direction  will  be  reselected,  if  the  selected  ray 
intercepts  the  sphere  of  influence,  or  "critical  sphere",  around  a "live"  de- 
tector, i.e.,  a detector  whose  scoring  capacity  is  currently  active.  This  an- 
gular reselection  procedure  comprises  three  distinct  stages:  (1)  direction 

checking;  (2)  reselection  of  angle:  (3)  weight  adjustment. 


Tt  .i  "FUJP"  latent  is  PICK-ed  (J1234S=5) , FLUP  is  called  instead  of  CARP  . 


In  the  direction  checking  stage,  the  orientation  of  live  detectors  with 


respect  to  the  originally  selected  ray  is  analyzed.  If  the  ray  does  not  in- 
tersect any  detectors,  reselection  is  bypassed.  If  an  intersection  is  found, 
reselection  may  be  necessary,  but  prior  to  reselection,  potential  "conflicts" 
must  be  resolved.  Two  types  of  conflicts  are  possible:  (1)  the  intersected 

critical  sphere  overlaps  the  critical  sphere  of  another  live  detector; 

(2)  the  cone,  with  vertex  at  the  collision  (or  source)  point,  which  is  tangent 
to  the  intersected  critical  sphere,  overlaps  the  corresponding  cone  of  another 
live  detector.  A conflict  is  resolved  by  Russian  roulette,  in  which  all  de- 
tectors, but  one,  are  deactivated. 

Subsequent  to  the  direction  checking  stage,  at  most  one  live  detector 
sphere  is  intersected  by  the  originally  selected  ray.  If  such  an  intersection 
has  survived,  a new  direction  is  selected  in  the  reselection  stage,  as  described 
in  Appendix  E. 

The  final  stage  in  the  reselection  procedure  involves  the  calculation  of  a 
w i iiustment  factor  which  compensates  for  the  biasing  introduced  by  re- 

el -action 

. ui  oui.iin.-  CARLO 

is  th  Monte  Carlo  calculation  routine  It  controls  importance 
, , collision  mechanics,  and  scoring. 

■nil  source  particle  is  accompanied  by  its  associated  array  of  problem 
weights,  F , i= 1 , 2 , . . . , NPROB,  where  NPROB  is  the  number  of  problems. 

Source  particles  generated  by  subroutine  SOUPIC,  as  well  as  latents,  are 
tested  on  energy.  If  the  energy  is  below  the  lower  bound  of  the  current  super- 
tin  it  is  stored  as  a latent.  Otherwise  it  is  transmitted  to  CARLO,  which 
tra  k.s,  distributes  further  collisions,  if  any,  and  scores  answers 
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When  point  detectors  are  being  used,  and  if  the  extended  particle  path 


passes  sufficiently  close  to  a specified  detector,  preliminary  geometrical 
calculations  necessary  to  bias  the  collision  positions  are  carried  out  This 
biasing  is  required  to  make  flux  estimates  for  point  detectors  bounded 

Given  the  region  number  IR,  the  energy  E,  source  particle  position  vector 
XB,  and  the  direction  of  flight  WB,  a sampling  weight  w=wIR-wE-wwB  is  calculated. 

The  subroutine  DR1  is  called,  which  provides  the  total  macroscopic  cross 
sections  ip,  i=l , . . ,NPROB,  for  each  of  the  problems. 

The  geometry  routine  GltS^)  is  called  to  provide  the  distance  to  tht 
first  region  boundary  encountered  from  XB  in  the  direction  WB.  G1  also  pro- 
duces IR',  the  region  number  on  the  other  side  of  the  boundary,  and  X’,  the 
point  of  intersection  of  the  track  with  the  boundary.  If  IR  is  a scoring  re- 
gion, the  contribution  0^  to  the  flux  is  calculated  and  stored  for  all  problems, 
where 

0.  = (F.  -F.e“MiSl)— , i=l , . . . , NPROB. 
ill  p ^ 

Once  the  tracking  and  scoring  (if  called  for)  are  completed,  the  Monte 
Carlo  procedure  necessary  to  generate  collision  positions  is  invoked.  This 
procedure  uses  three  different  algorithms,  depending  on  whether  or  not  the 
special  biasing  for  bounded  estimation  of  point  detectors  is  needed,  where  the 
special  biasing  can  use  either  of  two  algorithms 

The  general  principle  used  for  the  selection  of  collision  positions  in- 
volves an  implicit  splitting  and  Russian  roulette  procedure.  Specifically, 
to  select  the  positions,  an  uniio rmal ized  probability  density  f (s)  is  defined, 
where  ; is  distance  along  the  track,  and  the  expected  total  number  of  collisions 
is  given  by  r f (s)ds  The  actual  number  of  collisions  is  then  either  |I] 

or  [I)+l,  where  the  probability  of  the  additional  collisions  is  I-|I]r  where 
(xl  the  greatest  integer  less  than  or  equal  to  x. 
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In  carrying  out  the  collision  dropping  procedure,  f(s)  is  treated  in 
a piecewise  fashion,  where  in  each  interval  f(s)  and  its  indefinite  integral 
have  the  desired  analytical  properties,  most  important  being  that  the  in- 
definite integral  be  easily  invertible.  Finally,  to  compensate  for  any 
biasing  involved  in  using  f (s) , the  particle  weights  {F^}  must  be  adjusted 
accordingly. 

The  remainder  of  the  procedure  involves  the  definition  of  f(s),  which 
may  assume  three  forms,  as  a consequence  of  bounded  biasing.  The  algorithms 
which  have  been  implemented  are  described  in  Appendix  E. 

If  a collision  is  dropped  by  means  of  the  above-described  procedure,  sub- 
routine CARSCA  is  called  to  complete  the  remainder  of  the  processing.  When 
all  the  collision  positions  have  been  selected,  the  subroutine  proceeds,  as 
follows. 

A test  is  made  to  determine  whether  the  next  region  is  a transmission 

* 

region  (if  it  is,  the  coordinates  S',  IR' , energy,  time,  etc.,  are  written 

on  > , and  whether  it  is  the  escape  region.  If  it  is  not  the  escape  region, 

particle  is  moved  to  the  boundary  by  setting  IR=IR' , by  computing  the  new 

W1 

legion  weight  W , and  by  setting  F.=(— — ),  where  W is  the  old  region  weight. 

I l W0  1 

We  again  define  FMAX  = the  largest  of  F 's  and  play  another  game  of 

If  FMAX<F  , the  particle  is  either  killed,  or  it  survives  with  the 
— z 

' . renormalized  by  dividing  all  of  them  by  FMAX.  If  FMAX>F  , DR1  is  called 
to  provide  the  new  total  macroscopic  cross  sections  for  all  problems  and  con- 
troJ  is  transferred  to  the  section  which  calls  the  geometry  routine  G1 (S  ) . 
the  tracking  continues  until  a kill  occurs  or  the  escape  region  is  reached. 

* 

Not  used  with  point  detectors. 
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After  the  tracking  is  completed,  the  latent  table  of  particles  within 


the  current  energy  superbin  is  examined.  If  the  table  is  not  empty,  the  last 
latent  is  picked  up  and  processing  continues.  (Latents  are  created  in  CARSCAi 
The  subroutine  is  terminated  when  the  latent  table  is  exhausted. 

Subroutine  CARSCA 

This  routine  controls  all  the  processing  necessary  after  a collision 
position  is  selected  in  CARLO  The  particle  weights  are  adjusted  to  compen- 
sate for  the  biasing  used  in  the  sampling.  Subroutine  DR3  is  called  to  carry 
out  the  scattering  mechanics,  producing  new  energy  and  direction. 

If  the  scattering  event  was  absorption,  or  if  the  energy  was  below  the 
cutoff,  the  procedure  is  terminated  If  not,  a sampling  weight  is  calculated 
end  Russian  roulette  is  played  if  the  comparison  with  FZ  warrants  it.  Finally, 
if  not  terminated,  the  particle  is  stored  as  a latent  in  one  of  two  lists,  de- 
pending on  whether  or  not  the  energy  is  within  the  energy  superbin  being  con- 
sidered. Moreover,  if  it  falls  within  the  current  superbin,  a Russian  roulettt 
procedure  may  be  necessary,  if  the  list  is  full,  in  order  to  make  space  avail- 
able The  procedure  terminates  after  completion  of  latent  storing. 

If  an  interaction  tape  is  being  created,  for  subsequent  use  as  a sour 
of  secondary  radiation,  CARSCA  will  oversee  the  proper  storing  of  the  para- 
meters of  the  interaction  event.  This  includes  both  absorption  and  scatter  in-  > 
events . 

Subroutine  DR1 

This  routine  computes  the  unperturbed  and  perturbed  macroscopic  total 
ro;s  sections  at  a given  energy  for  a given  composition  It  first  picks  uj 

■ill  the  relevant  perturbations  (IP's  of  type  1,  2 and  3 of  relevant  elements) 
it  'ii  DRl  call  subroutine  ASSIGN  to  assign  the  perturbations  to  each  problem 
■ i i ■ <d:  to  compute  the  macroscopic  total  cross  section  for  each  problem 


according  to  the  perturbation  description  of  each  problem.  When  called 


by  CARLO,  DRl  also  computes  the  sampling  probability  distribution  for  picking 
he  element  of  reaction.  The  latter  computations  are  bypassed  when  DRl  is 
called  from  FLUP  or  TRALA  during  BFAP  estimation,  as  signaled  by  the  integer 
flag,  IFA=1. 

Subroutine  DR3 

Subroutine  DR3  is  called  by  CARSCA  to  perform  the  actual  collision 
mechanics.  It  returns  to  CARLO  with  a new  energy  (EPRIM) , a new  direction 
(WP) , the  cosine  of  the  scattering  angle  (CSTHT) , and  an  integer  (NCDB)  de- 
noting the  type  of  interaction  event.  The  types  of  interaction  are  listed 
below: 


NCDB  Interaction 

1 discrete  isotropic  inelastic  scattering 

2 discrete  anisotropic  inelastic  scattering 

3 isotropic  elastic  scattering 
anisotropic  elastic  scattering 

5 Klein-Nishina  gamma-ray  scattering 

6 absorption 

continuum  isotropic  inelastic  scattering 
8 continuum  anisotropic  inelastic  scattering 

DR3  first  picks  the  element  of  reaction  from  the  sampling  distribution 
imputed  in  DRl,  and  adjusts  the  perturbation  weight  for  each  problem.  It 
then  computes  the  sampling  distribution  for  the  selection  of  a reaction  type 
for  the  chosen  element. 

Subsequent  to  the  selection  of  a reaction  type,  control  is  passed  to 
ither  DR31  (for  elastic  scattering)  or  DR32  (for  inelastic  scattering)  for 
• hi  selection  of  post-scattering  direction  and  energy. 


Subroutine  DR31 


This  routine  is  called  by  DR3  subsequent  to  the  selection  of  an  elastic 
scattering  event.  Its  basic  function  is  to  compute  a post-scattering  direc- 
tion and  energy,  In  addition,  if  point  detectors  have  been  specified 
input,  the  BFAP  estimation  procedure  is  initiated  by  calling  FLUF. 

When  a perturbation  of  angular  distribution  exists  (type  2 or  6)  at 
the  collision  energy  E,  DR31  will  call  SAMPIC  to  pick  a value  of  CHI  fror,  Uu 
sampling  CHI-table  (if  it  exists).  SAMPIC  adjusts  the  problem-dependent 
turbation  weights  automatically.  If  a sampling  table  is  not  available,  i ; 
will  pick  a CHI  from  the  unperturbed  distribution,  and  call  XWADJ  to  ad , 
the  perturbation  weights. 

Finally,  if  point  detectors  have  been  specified,  the  angular  reselectici 
procedure  is  initiated  with  a call  to  DR34. 

Subroutine  DR32 

This  routine  is  called  by  DR3  if  an  inelastic  scattering  event  hac 
selected.  Its  basic  functions  are  the  same  as  those  of  DR31  for  elastic- 
scattering,  including  FLUP  and  DR34  calls  for  BFAP  estimation  and  angle  i 
selection,  in  the  event  that  point  detectors  have  been  specified. 

When  a perturbation  of  secondary  energy  distribution  exists  (type  . 
at  the  collision  energy  E,  DR32  will  call  SAMPIE  to  pick  a value  of  ENN 
the  sampling  ENN-table  (if  it  exists).  SAMPIE  adjusts  the  problem- depone 
perturbation  weights  automatically.  If  a sampling  table  is  not  avail  a:  l< 
selection  is  made  on  the  basis  of  the  unperturbed  data,  and  perturbation 
weight  adjustment  is  performed  within  this  routine. 
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Subroutine  FLUP 


This  routine  performs  the  BFAP  estimation  for  each  live  detector.  It  is 
called  by  AMCAR  for  each  source  particle  (the  primary  source  must  be  isotropic) 
and  for  "FLUP  latents"  (i.e.,  pseudo  particles,  established  in  a previous  call 
to  FLUP  for  estimation  tracking,  whose  energy  fell  below  the  then  current  super- 
bin cutoff).  It  is  also  called  by  DR3,  DR31,  and  DR32  for  Klein-Nishina,  elastic, 
and  inelastic  collisions,  respectively. 

The  estimator  scored  is  an  array,  each  element  of  which  is  given  by 

-S. 

W^g  (cosQ ) e 1/4irr  , i=l,...,  n 
where  n = the  number  of  problems; 

W = net  particle  weight,  adjusted  for  all  sampling  biasing  that 
may  have  taken  place  previously; 

g (cosG)  = probability  of  scattering  through  an  angle  of  cosine 
9,  (or  1 for  an  isotropic  source); 

S = optical  distance  between  collision  (or  source)  point  and 
etector  point; 

r = problem  independent,  geometric  equivalent  of  S. 

At  a collision  point,  the  g^array  is  computed  by  calling  DR33,  which 
i returns  an  energy  EP,  the  lab  energy  resulting  from  a pseudo -scattering 
the  direction  of  the  detector.  If  EP  is  below  the  current  superbin  cutoff, 
i i ! the*  information  already  computed  for  the  estimator  ray  is  stored  (via  calls 

0 PACK  and  PICK),  to  be  picked-up  later  as  a "FLUP  latent".  The  S . -array  is 
computed  via  subroutine  TRALA. 

The  detectors,  for  which  a flux  estimator  is  computed,  are  determined 

1 fiLlows;  For  a "FLUP  latent",  only  the  detector  which  was  being  treated  at 
ae  time  of  the  latent  store  is  treated.  Otherwise,  the  variable  KDLIV  (see 

vppendix  F)  is  examined.  For  KDLIV^O,  all  detectors  are  "alive",  i.e.,  their 


scoring  capacity  is  active.  But  for  KDLIV<0,  only  detector  # -KDLIV  is 


"alive",  as  a result  of  Russian  roulette,  which  was  performed  to  resolve  a 
reselection  "conflict"  (see  Appendix  F and  the  discussion  of  Subroutine 
ARG) . In  the  latter  case,  only  detector  # -KDLIV  is  treated,  with  an  addi- 
tional weight  adjustment  factor  NDEf  (the  total  number  of  detectors),  which 
compensates  for  the  Russian  roulette  biasing. 

Subroutine  PACK  (X,WX ,E , IR,T , IDET,W,NHIST,WC, J12345 ,F ,P) 

This  subroutine  transfers  the  items  in  the  calling  sequence  to  the  P 
array.  The  relationship  between  the  arguments  and  the  array  P is  shown  t • 1 . 


Argument 

P 

Definition 

X (1) 

P(l) 

X (2 ) 

P(2) 

Position  Coordinates 

X (3) 

P (3) 

WX(1) 

P(4) 

WX  (2) 

P(5) 

Direction  Cosines 

WX  (3) 

P<6) 

E 

P(7) 

Energy 

IR 

P(8) 

Region  Number 

T 

P(9) 

Time 

IDET 

P(10) 

Detector  Number 

W 

P(ll) 

Sampling  Weight  = 

^IR'V'W 

NHIST 

P (12) 

History  Number 

WC 

P (13) 

Weight  (carry-along 

weight,  usually  =1.) 

.712345 

P (14) 

Particle  Type  (defined 

below  in  description 

in  PICK  routine) 

F. 

l 

P(14+i) 

Problem-dependent  weights. 

i=l 10 
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Subroutine  PICK 


We  have  seen  throughout  the  previous  sections  that  particles  degraded 
below  the  energy  EBL,  low-energy  limit  of  the  superbin  currently  treated, 
were  stored  as  latents  by  calling  the  subroutine  PICK.  They  were  later  picked 
up  by  AMCAR  by  calling  the  subroutine  PICK.  The  subroutine  INPUTD  allocates 
the  memory  to  data,  scores,  etc.  The  remaining  memory  is  assigned  to  the  sub- 
routine PICK,  to  be  used  as  a buffer  for  latents.  One  'end'  of  the  buffer  is 
assigned  to  'degraded'  particles.  This  is  the  end  of  the  buffer  where  particles 
are  being  stored.  The  other  'end'  of  the  buffer  is  assigned  to  'unsorted' 
particles,  i.e.,  the  particles  to  be  picked.  Associated  with  each  end  of  the 
buffer  is  a magnetic  tape  to  be  used  when  the  buffer  overflows.  There  are  two 
modes  of  operation.  In  one  mode,  the  top  of  the  buffer  is  unsorted  and  the 
bottom  is  sorted.  When  the  switch  is  made  from  one  superbin  to  the  next,  the 
'unsorted'  part  is  empty,  and  the  'degraded'  part  may  have  particles  which 
become  'unsorted'  for  the  supergroup  about  to  be  treated.  The  designation  of 
the  buffe's  (and  of  the  tapes)  is,  therefore,  switched. 

There  is  no  set  boundary  between  the  two  ’ends'  of  the  buffers.  The 
number  of  particles  in  the  'unsorted'  buffer  keeps  decreasing,  whereas  the 
number  of  particles  in  the  'degraded'  buffer  keeps  increasing,  and  can  increase 
taster  than  the  other  number  decreases.  Therefore,  the  two  parts  of  the 
buffer  can  meet,  causing  an  overflow  of  the  buffer. 

It  is  then  determined  which  'end'  of  the  buffer  is  longest,  and  a number 
of  particles  exactly  equal  to  one-half  the  total  length  of  the  buffer  are 
written  from  the  longest  'end'  onto  the  corresponding  magnetic  tape. 
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When  the  'unsorted'  buffer  becomes  empty,  a test  is  made  to  determine 
if  any  'unsorted'  particles  are  available  on  the  corresponding  magnetic  tape. 

If  none  are  available,  the  calculation  has  been  completed  for  the  current 
superbin.  If  some  are  available,  they  are  read  into  the  buffer  if  room  is 
available.  If  room  is  not  available,  it  is  made  available  by  writing  out  part 
of  the  other  buffer  on  the  other  tape;  the  length  of  the  record  written  out 
from  one  'end'  is  equal  to  the  length  of  the  record  to  be  read  into  the 
other  ' end' . 

The  subroutine  PICK  deals  with  different  kinds  of  latents.  The  quantities 
stored  are;  X,fi,E,IR,T,I,W,NHIST,WC,F^’s  and  J12345,  where  X is  the  position, 

£2  the  direction,  E the  energy,  IR  the  region  number,  T the  time,  W and  F^'s 
the  weights,  NHIST  the  history  number,  and  WC  a normalization  factor.  I and 
J12345  are  indices. 

J12345  = 1 identifies  a source  particle 

= 2 identifies  a particle  coming  out  of 
elastic  scattering. 

= 3 identifies  a particle  coming  out  of 
inelastic  scattering, 

= 5 identifies  a FLUP  latent 

(In  other  parts  of  the  code,  J12345=6  identifies  a transmitted  particle, 
J12345=7  an  inelastic  interaction,  J12345=8  an  absorption,  J12345=9  an  elastic 
interaction,  and  J12345=10  a non-elastic  interaction.) 

Subroutine  SEEK  (E, EOUT, NOUT, I) 

Given  the  vector  array  EOUT,  of  length  NOUT+1,  and  the  argument  'E'  the 
routine  will  perform  a binary  search  and  return  'I'  such  that 


EOUT ( I ) <E<EOUT ( 1+1 ) . 


Subroutine  SOUCAL 


! 


SOUCAL  will  read  source  data  and  prepare  tables  for  use  by  SOUPIC. 

The  first  card  contains: 

NSR  Number  of  regions  where  the  source  extends.  If  this  number 

is  0,  an  external  source  tape  is  expected,  and  no  further 
input  is  required. 

IFLAG  Number  of  energies  used  to  define  the  source  spectrum.  The 
spectrum  is  specified  later  in  input  which  is  in  the  form 
of  a histogram.  It  will  be  defined  by  IFLAG  entries. 

ISW  Switch  determining  the  normalization  of  problem.  If  it  is 

0,  fluxes  will  be  normalized  to  a single  (unbiased)  source 
particle.  If  it  is  1,  fluxes  will  be  normalized  to  the 
'total  power',  i.e.,  output  will  be  flux  per  source  particle, 
multiplied  by  the  summation  (over  source  regions)  of  power 
density  times  volume. 

This  is  followed  by  (NSR)  cards  specifying  the  (NSR)  source  regions. 

Each  of  these  cards  gives: 

ISR  Specification  of  a geometrical  region  number. 

P Power  density  in  that  region. 

ISO  ISO=0,  isotropic, 

ISO=l,  monodirectional.  If  ISO=l  in  any  of  the  (NSR) 
regions,  a monodirectional  source  will  be  assumed  in  all 
the  (NSR)  regions. 

After  this  follow  the  (IFLAG)  cards  specifying  the  energy  spectrum  assumed 
to  apply  to  all  the  source  regions.  This  is  a table  of  E vs.  S(E)  «*  the  height 


of  the  histogram  to  the  right  of  the  energy  E.  SOUCAL  will  compute  the  integral 
spectrum 

F (E) 
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Linear  interpolation  is  assumed  on  E vs.  F(E).  The  first  entry  must  be  for 

E<E  .,  and  the  last  for  E>E  . . 

— cut  — high 

If  there  is  time  dependence  (NT>0),  the  time  dependence  of  the  source  must 
be  specified.  A card  gives: 

NOT  Number  of  cards.  This  card  is  followed  by  NOT 

cards  giving 

r 

t vs.  I S^ (t) dt. 


Linear  interpolation  is  assumed  between  entries 
in  the  table. 

Finally,  if  the  source  is  monodirectional,  the  projections  of  the 

direction  must  be  given  on  the  last  card. 

The  subroutine  SOUCAL  reads  in  all  this  input,  prints  it  back,  and  pre-com- 
putes tables  to  pick  directly  from  the  biased  source  distribution.  Since  we 
allow  perturbations  of  the  source  spectrum  which  are  input  as  density  func- 
tions in  the  form  of  histograms,  two  preliminary  steps  are  necessary  before 
have  the  final  biased  source  distribution  for  subroutine  SOUPIC: 

1.  Form  the  envelope  of  the  unperturbed  and  perturbed  spectra, 

2.  Calculate  the  distribution  function  of  this  envelope. 

Step  (1) 


a.  Normalize  each  source  spectrum  between  E . and  E,  . , . 

cut  high 

b„  Construct  a table  containing  all  the  energy  boundaries 
present  in  the  specifications  of  the  various  source 
spectra  (histograms) . 

c.  Rearrange  this  table  into  a monotonic  increasing  sequence. 

d.  Remove  duplicate  energies  from  the  ordered  sequence.  This 
provides  a combined  energy  table  for  calculating  the 
envelope  spectrum. 
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e.  For  each  energy  E surviving  step  d,  call  SEEK  to  locate  it  in 
each  of  the  input  spectrum  energy-tables,  and  calculate  the 
strength  of  each  source  spectrum  (dp/dE)  at  this  energy. 

f.  For  each  of  these  energies  E,  take  the  maximum  individual 
spectrum  strength  as  the  envelope  value  corresponding  to  the 
energy  range  from  E to  the  next  higher  E in  the  combined 
energy  table.  Since  the  individual  spectra  are  assumed  to 
have  piecewise  constant  dp/dE,  and  the  combined  energy  table 
includes  all  of  the  histogram  energies,  no  details  of  any 
spectrum-structure  are  lost  in  this  construction  of  the  en- 
velope. 

Step  (2) 

Calculate  the  envelope  integral  spectrum  and  normalize  it,  and  we  have  the 
sampling  source  spectrum  S (E) . 

The  code  first  pre-computes  a table  of 

r - 

SPEC (I)  = / S (E)dE, 

7ei 


where  the  E^'s  take  the  values  of  of  all  the  energy  boundaries  where  the 

energy  weight  changes,  and  Ecut»  The  code  then  runs  through  all  the  source  re- 
gions and,  for  each  new  energy-importance  set  encountered,  pre-computes  a table  of 


SPEC (I, J)  = 


f 


J 


S (E)dE 

VE) 


where  J runs  from  1 to  the  total  number  of  different  energy- importance  sets  en- 
countered in  the  source  regions;  also  for  each  new  angular  importance  set  en- 
countered, a table  f 


P< 


/ 1 

I,K)  = / 


dto 


I W (u>) 

I U).  '0 
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where  the  w^'s  take  the  values  of  cose  at  which  the  angular  weight  changes, 
and  K runs  from  1 to  the  total  number  of  different  angular  importance  sets 
encountered.  The  different  tables  are  renormalized  and  both  the  modified 
and  unmodified  integrated  source  are  computed  in  each  source  region.  The 
former  quantities  are  proportional  to  the  probability  with  which  particles 
should  be  picked  in  different  source  regions.  A table  SOUR(L)  is  built  up, 
which  gives  the  cumulative  probability  for  a source  to  be  picked  in  the  Mth 
region  for  M£L. 

Subroutine  SOUPIC 

This  is  u subroutine  which  picks  particles  from  the  biased  source  dis- 
tribution. 

If  an  external  source  tape  is  to  be  used,  groups  of  35  source  particles 
are  read  from  tape  into  a buffer,  and  returned  one  by  one  to  the  main  code. 
The  quantities  describing  a source  particle  are: 


XB  Coordinates  of  the  particle 

IR  Region  number  where  particle  is  bom 

WB  Direction  of  the  particle 

T Time  at  which  particle  is  born 

E Energy  of  the  particle 

NHIST  History  number  attached  to  the  source  particle 

i=l,...,  total  number  of  problems.  Statistical 
weight  of  the  particle  (usually  set  equal  to 
unity  if  no  perturbation  of  the  source  spectrum) . 


The  procedure  for  internal  source  generation  is  outlined  below. 

A first  random  number  £ is  compared  to  the  table  SOUR(L)  (computed  by 
SOUCAL) . The  smallest  L for  which  SOUR(L)^  determines  the  region  IR=ISR(L) 
to  be  picked  from.  Standard  techniques  are  used  to  pick  coordinates  of  points 
uniformly  distributed  in  a region. 
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The  energy  is  picked  next.  A stratified  random  number  (called  CE 
in  the  code)  is  obtained  (stratification  is  done  for  each  statistical  aggre- 
gate of  source  particles) . A biased  random  number  £ is  then  obtained  by 
interpolation  in  the  SPEC  vs.  SPEC(J)  tables  pre-computed  by  SOUCAL.  (The 
J is  determined  by  the  region  number.)  Finally,  the  energy  is  determined  by 
solving  the  equation 


S (E)dE 


using  linear  interpolation. 


Figure  2 - Selection  of  £ for  Picking  a Source  Energy 

The  direction  WB  of  the  source  particle  is  determined  as  follows.  If 
the  source  is  isotropic  and  there  is  angular  importance  sampling,  the  cosine  of 
the  angle  between  the  particular  aiming  angle  and  the  direction  is  chosen  by 
picking  a random  number,  and  interpolating  between  the  angular  mesh  supplied 
on  input  vs  the  table  P(I,K)  pre-computed  by  SOUCAL.  (The  K is  determined 

by  the  region  number  ) A random  azimuth  is  then  picked,  which  completes  the 

* 

specif ication  of  the  direction  . In  the  absence  of  angular  importance  in  the 

* If  point  detectors  exist,  the  angular  reselection  procedure  is  initiated 
by  calling  ARG. 
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source  region,  standard  techniques  are  used.  The  case  of  a monodirectional 


source  also  can  be  handled,  provided  there  is  no  angular  importance  in  the 
source  region.  Finally,  if  time  dependence  is  to  be  determined,  a time  T 
is  determined  by  another  random  number  fi  and  the  solution  of  the  equation 


C = 


S (t) dt. 


The  quantities  communicated  to  the  main  code  are  XB,  IR,  WB,  T,  E,  NHIST,  W, 
and  where  (W^  = 1) , and  W is  the  weight  in  the  region  IR  at  energy  E in  the 
direction  WB. 


Subroutine  TALLY  (J, NHIST, NCOL,NDEG, NABS, FKILL, BIRTH, ESCAP,NRMAX) 

The  routine  is  used  to  print  a tally  at  the  end  of  each  statistical  aggre- 
gate and  to  print  a tally  by  region  at  the  end  of  the  problem. 

At  the  end  of  each  aggregate  the  following  items  are  printed. 

NC  Total  number  of  collisions  thus  far 

ND  Total  number  of  degrades  thus  far 

NA  Total  number  of  absorptions  thus  far 

FK  Total  number  of  kills  thus  far 

BI  Total  number  of  births  thus  far 

ES  Total  number  of  escapes  thus  far 

T Elapsed  time  thus  far 

JCOUNT  Total  number  of  particles  on  the  interaction/ 
transmission  tape. 

The  same  items,  except  for  ES  and  T,  are  printed  at  the  end  of  the 
problem  as  a function  of  region. 
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Subroutine  TRALA 


This  routine  is  called  by  FLUP,  in  the  BFAP  estimation  procedure,  to 
compute  the  problem-dependent  array  of  optical  distances  from  collision  (or 
source)  point  to  a detector  point. 

Given  an  intial  point  XB  and  a direction  WB,  the  routine  will  track 
from  XB  in  the  direction  WB.  The  tracking  will  continue  until  the  geometric 
distance  to  the  detector  has  been  reached,  or,  as  a result  of  Russian  roulette, 
the  estimator  is  set  to  zero  to  reduce  unproductive  tracking  (see  Appendix  G) . 
Subroutine  UNPACK  (X ,W, E , IR,T, IDET, F ,NHIST, WC , J12345 , P ) 

Subroutine  UNPACK  distributes  the  24-word  P array  among  the  variables  in 
the  argument  list  according  to  the  correspondence  shown  in  the  discussion  of 
subroutine  PACK. 

5.4,2  Glossary  of  Important  Variable  Names 

The  following  is  a brief  outline  of  major  common  blocks  in  the  program. 

A detailed  description  of  each  variable  in  common  is  given  on  the  following 
pages. 

The  "blank"  common  block  contains  the  master  storage  array  and  frequently 
used  input  arrays. 

Common  REGPAR  contains  six  parameters  used  in  the  Monte  Carlo  calculations. 
Common  INPUT  contains  non-subscripted  input  parameters. 

Common  PAREM  contains  position,  direction,  energy  and  importance  sampling 

parameters. 

Common  COMPUT  contains  parameters  computed  from  input. 

Common  CROSA  contains  cross  section  parameters. 

Common  FAP  contains  f lux-at-a-point  parameters. 

Common  SOU  contains  source-dependent  parameters. 

Common  PUTAD  contains  data  controlling  the  parameters  to  be  written  on 
the  interaction,  transmission  tape. 
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Blank  Common 


EOUT(IOO) 


EWTAB (50) 
ANGLE (50) 
TTAB (50) 
ASTER (5000) 


An  array  containing  the  output  energy  bins  for 
flux  results.  The  array  contains  all  bins  for 
all  output  supergroups. 

The  energy  mesh  for  energy  importance  sampling. 
The  cosine  mesh  for  angular  importance  sampling. 
The  time  mesh  for  time -dependent  problems. 

The  master  storage  array  containing  input  and 
flux  data.  A complete  description  appears  in 
Appendix  C.  The  array  is  equivalenced  with  array 
MASTER. 


COMMON  REGPAR 

ISC  The  location,  in  MASTER,  of  the  flux  scores  for 

scoring  region  ISC. 

NREG  The  composition  number  for  region  IR. 

IRW  The  location,  in  MASTER,  of  the  region  weight 

for  region  IR. 

IEW  The  location,  in  MASTER,  of  the  energy  weight  table 

for  region  IR. 

I AM  The  location,  in  MASTER,  of  the  aiming  angle  for 

region  IR. 

IANG  The  location,  in  MASTER,  of  the  angular  weight 

table  for  region  IR. 

Note  that  the  above  six  items  are  all  keyed  to  a region  number,  "IR". 

The  six  items  are  retrieved  from  the  MASTER  array  by  a call  to  UNPR. 
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Common  INPUT 


NSTART 

NSTOP 

NSTAT 

NRMAX 

NG 

NT 

NOUT 

NUMSC 

NRWL 

I REX 

NEWL 

NEW 


NAIML 

NUMANL 

NUMANG 

JRT 

ECUT 

ETHERM 


The  number  of  real  time  seconds  of  running 
time  before  terminating  and  editing. 

Number  of  the  last  history  to  be  treated. 

Number  of  histories  per  statistical  group. 

Number  of  regions  in  the  geometry. 

Either  0 for  a neutron  problem  or  1 for  a 
gamma  problem. 

Number  of  output  time  bins. 

Number  of  output  energy  bins. 

Number  of  flux  scoring  regions. 

Number  of  distinct  region  weights. 

The  escape  region  number. 

Number  of  energy  bins  for  energy  weighting. 
Number  of  distinct  energy  weight  sets.  If 
NEWL  and  NEW  = 0,  the  problem  contains  no 
energy  weighting. 

Number  of  distinct  aiming  angles. 

Number  of  angular  bins  for  angular  weighting. 
Number  of  distinct  angular  weight  sets. 

Not  used. 

Low  energy  cutoff  (ev) . Tracking  of  a particle 
is  terminated  if  its  energy  degrades  below  ECUT. 
Thermal  energy  if  a thermal  group  is  required. 
ETHERM  must  be  within  the  energy  limits  of  the 
problem. 
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TCUT 


FZ 

EHIGH 


EBL 

EBH 

XB(3) 

WB  (3) 

E 

IR 

T 

IDET 

F 

NHXST 

WC 

J12345 

WP  (3) 

XP(3) 

EPRIM 

ATWT 

NCDB 

CSTHT 


Time  cutoff  (equal  to  first  entry  of  time 
bin  boundary  table). 

See  discussion  in  Section  5.4.1. 

High  energy  cutoff  (ev) . This  should  be  less 
than  or  equal  to  the  highest  energy  for  which 
cross  sections  are  available. 

Lower  bound  of  the  current  superbin. 

Upper  bound  of  the  current  superbin. 

Common  PAREM 

The  X,Y,Z  coordinates  of  the  current  particle's 
starting  point. 

The  direction  cosines  of  the  current  particle. 
Energy  of  the  particle. 

Region  number  of  the  particle. 

Time  of  flight  of  the  particle. 

A detector  number  used  in  f lux-at-a-point 
calculations. 

Importance  sampling  parameter. 

Current  history  number. 

A weight  parameter  calculated  by  SOUPIC  . 

A particle  type  flag  (see  the  discussion  of  PICK; 
Direction  cosines  of  particle  after  scatter. 

X,Y,Z  coordinates  of  current  position. 

Energy  after  scatter. 

Atomic  weight  of  scattering  element. 

Interation  type  indicator  (see  discussion  of  DR3) 
Cosine  of  scattering  angle. 


u 

LCHI 

XATWT 

IERR 

IDBG 

IRPRIM 

NASC 

LSURF 

XXX 

LRI 

LRO 

RIN 

ROUT 

KLOOP 

LOOP 

LTYPE 

PINF 

NOA 

DIST 

NUMNOU 


JONUM 

LNCNOL 


Total  macroscopic  cross  section  at  energy  E 
for  region  IR. 

Location  in  the  XS  array  of  a table  used 
in  anisotropic  scattering. 

Identification  digit  of  the  scattering  element. 
Error  indicator. 

Debug  printout  flag. 

Next  region  to  be  entered  by  the  ray. 

A flag  to  initiate  the  Gl  routine  for  a new  ray. 
Not  used. 

Not  used. 


Geometry  subroutine  parameters 


Common  COMPUT 

The  product  of  NUMSC  (the  number  of  scoring 
regions)  and  NOUT  (the  number  of  output  energy 
bins) . 

An  index  used  in  flux  scoring. 

The  location  in  MASTER  of  the  collision  by 
region  table. 


108 


LBIRTH 

LREGT 

LFKILL 

LESCAP 

LLAST 

NDQ 

LNDEG 

LNABS 

LSCORE 

LPACK 

NTOT 

LGEOM 

LEGEOM 

NB 

JUPPER 
J LOWER 
KPHYS 
NELEM 
LENCHI 


The  location  in  MASTER  of  the  birth  by  region 
table. 

The  location  in  MASTER  of  the  region  data  table. 
The  location  in  MASTER  of  the  kills  by  region 

table. 

The  location  in  MASTER  of  the  escapes  by  region 
table. 

The  location  of  last  word  in  MASTER  used  by 
the  program. 

The  size  of  the  MASTER  array  (set  in  SAMCAR) . 

The  location  in  MASTER  of  the  degrade  by  region 
table. 

The  location  in  MASTER  of  the  absorption  by 
region  table. 

The  location  in  MASTER  of  the  flux  scoring  array. 
Not  used. 

NUMSC  times  (the  number  of  energy  bins  in  the 
largest  supergroup). 

The  lo  tion  in  MASTER  of  geometry  data. 

The  last  location  in  MASTER  of  geometry  d^lu. 
Common  CROSA 

Cross  section  band  counter. 

Indices  in  the  ETOT  array  giving  the  energy 
range  of  the  current  supergroup. 

Number  of  physical  compositions. 

Number  of  unique  elements. 

Length  of  the  CHI— table  for  anisotropic 
scattering . 
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LENENN 


Length  of  the  ENN-table  for  inelastic 


NLOC 


NENERG 


NENBAN 
NBAND 
LENXS 
EBLX 
NAME (30) 
ATTAB (30) 
LOCS (30) 
LOCLEV (30) 

I.ENPLV  ,30) 
ETOT ( 500) 
XS (5000) 

XAD (25) 

YAD (25) 

ZAD (25) 

ID 

NDFAP 

NDET 

LSCFAP 

LPAFAP 


continuum  scattering. 

Number  of  words  in  the  XS  array  for  current 
supergroup. 

Total  number  of  energies  in  the  cross  section 
mesh. 

Number  of  energies  in  the  current  supergroup. 
Number  of  cross  section  bands. 

Size  of  the  XS  array 

Lower  energy  bound  of  the  bands. 

Array  of  element  identifiers. 

Atomic  weight  for  each  element. 

Location  in  XS  of  o for  each  element. 
Location  of  the  ELEV-table  for  inelastic 
discrete  scattering. 

Length  of  the  PLEV-taL) , . 

The  energy  mesho 

The  cross  section  data  array. 

Common  FAP 

The  array  of  X-coordinates  for  the  detectors. 
The  array  of  Y-coordinates  for  the  detectors. 
The  array  of  Z-coordinates  for  the  detectors. 
The  number  of  the  detector  being  processed. 
Not  used. 

The  total  number  of  detectors. 

Not  used 
Not  used. 
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Common  PUTAD 


t 


IRT (10) 

XWA 

IWI 

IWE 

IWD 

IWO 

OTHER (5) 


Transmission  region  numbers.  Only  the  first 
three  are  used. 

Not  used. 

Not  used. 

Elastic  recording  flag.  Nonzero  for  recording. 
Not  used. 

Non-elastic  recording  flag.  Nonzero  for 
recording . 

Not  used. 


5.4,3  Description  of  Output 

The  printed  output  consists  of  four  parts: 

1.  Title  card  and  cross  section  band  energy  limits, 

2.  Geometry  data, 

3.  Monte  Carlo  data, 

4.  Results  of  the  Monte  Carlo  calculation. 

These  four  parts  are  discussed  below. 

1,  Cross  Section  Data 

The  printed  lines  are  basically  a repeat  of  input  items  3 and  5 (see 
section  5.3). 


2.  Geometry  Data 

The  printed  geometry  output  is  simply  a repeat  of  input  item  7 (see 
section  5.3).  In  addition,  the  storage  requirements  in  the  geometry  arrays 
are  also  given. 

3.  Monte  Carlo  Data 

As  in  Geometry  Data,  the  Monte  Carlo  data  are  a repeat  of  the  input 
data  (items  8 through  30). 


Ill 


4.  Results  of  the  Monte  Carlo  Calculation 


The  program  will  print  for  each  statistical  aggregate  the  following 
i terns : 

History  number 

Number  of  collisions 

Number  of  degrades 

Number  of  Absorptions 

Number  of  kills 

Number  of  births 

Number  of  escapes 

The  elapsed  real  time  in  seconds. 

The  number  of  transmissions  and  interactions  recorded  on  tape. 

Note  that  all  of  the  above  items  are  cumulative  and  are  printed  as  a 
single  line  for  each  aggregate. 

When  all  the  histories  have  been  prcc’s-rd,  a printout  giving  t?  • to  ai 
number  of  absorptions,  elastics,  and  transmissions  will  occur. 

The  previously  discussed  tallies  as  a function  of  aggregate  will  be 
repeated  as  a function  of  region.  Thus,  for  each  region  the  number  of 
collisions,  degrades,  absorptions,  etc.,  will  be  printed. 

5.5  Tape  Utilization 

The  following  describes  the  function  of  each  tape  used  in  the  program. 
Tape  numbers  refer  to  FORTRAN  logical  numbers.  All  tapes  are  used  in  the 
binary  mode. 

Tape  8 

The  processed  perturbation  data  tape.  (See  Section  4.4). 

Tape  10 

The  organized  data  tape  (ODT) . The  tape  contains  cross  section  data 
for  a given  problem. 


112 


Tape  1 2 


Sampling  CHI-  and/or  ENN-tables,  generated  by  SAMSAM.  (See  Section  4.4). 
Tape  14 

The  interaction/transmission  tape.  All  interactions  and  transmission; 
are  written  on  this  tape  for  use  in  subsequent  problems. 

Tape  15 

An  external  source  tape.  The  tape  may  also  be  the  transmission  part, 
of  tape  14. 

Tape  16 

The  statistical  aggregate  tape.  The  AMCAR  routine  uses  this  tape  it 
record  each  completed  aggregate.  The  edit  routines  then  process  the  tape 
to  obtain  the  final  flux,  dose  results. 

Tape  17  & Tape  18 

Temporary  storage  tapes  for  latents.  The  tapes  are  used  by  the  PICK 
routine. 

* 

■>.6  Error  Messages 

Three  types  of  error  indications  are  given  by  the  SAMCEP  program. 

Type  1 errors  give  an  error  message  and  cause  the  program  to  terminate. 

Type  2 errors  give  an  error  message,  but  do  not  terminate  the  calculation. 

Type  3 errors  terminate  the  calculation,  but  give  no  error  message.  The 
possible  error  stops  and  messages  are  discussed  below. 


* 

A discussion  of  HAND  messages  is  included  for  completeness,. 


1 1 1 


Type  1 Errors 


1.  SEEK  ERROR 

The  error  occurs  in  SEEK  and  is  caused  by  an  argument  out  of 
range  of  the  vector  being  searched.  The  most  probable  causes 
are  an  input  error  in  the  output  energy  bins  or  the  energy  or 
angular  weight  bins. 

2.  OUT  OF  RANGE  ON  EBAND 

The  error  occurs  in  BAND  and  is  caused  by  the  input  cross 
section  energy  band  mesh  being  outside  the  range  of  the  energies 
on  the  EDT.  The  most  probable  causes  are  an  input  error  or  using 
the  wrong  EDT. 

3.  NO  MORE  ELEMENTS  ERROR  IN  BAND 

The  error  occurs  in  BAND  and  is  given  when  an  isotope  identifier 
in  the  BAND  input  cannot  be  found  on  the  EDT.  The  most  probable 
causes  are  an  input  error  or  using  the  wrong  EDT. 

4.  ***  ERROR  - BOTH  FIRST  AND  LAST  BIN  BOUNDARIES  MUST  BE 
FLAGGED  WITH  NEGATIVE  SIGNS 

5.  ***  ERROR  - EHIGH  MUST  BE  WITHIN  ENERGY  BINS 

6.  ****  ERROR  - ECUT  MUST  BE  WITHIN  ENERGY  BINS 

7.  ****  ERROR  - ET1IRM  MUST  BE  WITHIN  ENERGY  BINS 

8.  THE  NUMBER  OF  ENERGY  BINS  IS  TOO  BIG  THE  MAX  IS 

The  error  indicates  too  many  energy  output  bins  to  fit  in  the 
machine.  The  input  must  be  modified.  A suggestion  is  to  allow 
another  output  supergroup. 

9.  ERROR  IN  NUMSC 

The  error  occurs  when  a scoring  region,  as  supplied  by  the  region 
parameter  input,  is  bigger  than  the  input  value  of  NUMSC.  The 
input  should  be  checked. 
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10.  NO  ROOM  FOR  DATA 


The  error  indicates  that  the  total  room  occupied  by  the  input 
is  greater  than  the  allowable  maximum.  The  input  must  be 
modified. 

Note  that  errors  (4)  through  (10)  all  occur  in  INPUTD.  These  errors 
will  allow  the  remainder  of  the  input  to  be  processed  but  the  program  will 
not  perform  the  Monte  Carlo  calculation. 

11.  ****  ERROR  - SPECTRUM  NOT  DEFINED  BETWEEN  EHIGH  AND  ECUT 
The  error  occurs  in  SOUCAL  and  indicates  that  the  energy  spectrum 
violated  one  of  the  following  constraints: 

ECUT  >_  ET  (1 ) 

EHIGH  <_  ET(IFLAG) 

where  ET  is  the  input  energy  spectrum  and  IFLAG  is  the  number  of 
entries  in  the  spectrum  table. 

12.  CANNOT  HAVE  ANG.  IMP.  FOR  ANISOTROPIC  SOURCE 

The  error  occurs  in  SOUCAL  and  indicates  that  the  special  source 
direction  option  was  used  in  a region  containing  angular  importance. 
This  condition  is  not  allowed. 

13.  ERROR  IN  NOUT 

The  error  indicates  that  a source  energy  was  calculated  by  SOUPIC 
and  is  outside  the  range  of  the  output  energy  mesh.  Check  the  input. 

14.  , BAD  IRPRIM  IN  CARLO 

An  IRPRIM  of  zero  was  calculated  in  the  G1  routine.  The  geometry 
input  and  the  source  position  data  should  be  checked. 
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15.  SI  OUTSIDE  BOUNDS 


The  error  indicates  that  SI,  the  distance  to  the  next  boundary, 
as  calculated  by  Gl,  is  either  zero  or  greater  than  10^.  The 
geometry  input  should  be  checked.  A more  serious  cause  is  an 
error  in  compilation. 

16.  ERROR  IN  SP  IN  CARLO 

"S"  the  distance  to  the  next  collision  is  greater  than  SI, 
the  total  distance  through  the  region.  The  FORTRAN  statements 
in  CARLO  should  be  checked  for  a compile  error. 

17.  ERROR  IN  NREG 

A composition  number  of  zero  or  greater  than  the  KPHYS,  the 
number  of  compositions,  was  encountered.  Check  the  region  input. 

18.  ERROR  IN  IRPRIM 

A region  number,  IRPRIM,  was  greater  than  the  number  of  regions. 
Check  the  region  input. 

19.  ERROR  IN  ISC 

A scoring  region  number,  ISC,  was  greater  than  the  number  of 
scoring  regions.  Check  the  region  input. 

20.  ERROR  IN  NCDB  IN  CARLO 

An  illegal  interaction  digit,  NCDB,  has  been  generated  by  DR3, 
the  interaction  routine.  The  cross  section  data  should  be 

checked. 
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6.  PROGRAM  SAMOUT 


6.1  General  Discussion 
— 

Program  SAMOUT  reads  the  scores  for  each  statistical  aggregate  in  super- 
groups from  the  output  tape  (IAGG=tape  16)  of  the  previous  program,  SAMCAR,  and 
writes  a final  cumulative  statistical  tape  (NAGG=tape  4 in  this  program)  from 
which  all  editing  will  be  done.  Subroutine  EDIT  will  print  the  normalized  flux 
with  the  percentage  error  for  all  the  problems,  and  the  differences  (with  their 
percentage  errors)  between  all  or  any  two  problems.  Once  the  cumulative  statis- 
tical tape  (NAGG)  is  written,  program  SAMOUT  can  be  rerun  using  only  the  NAGG 
tape  to  edit  any  other  problems  or  differences  desired. 

Response  functions  (cf.  Section  5.2.7)  will  be  read  in  at  the  time  the 
IAGG  tape  is  processed;  the  response  functions  of  1.0  are  built  into  the  program. 

6.2  Description  of  the  Flux  and  Percentage  Error  Calculation 
The  output  of  SAMCAR  consists  of  flux  integrals 

F(IR,IE,IT,IA) 

* 

integrated  over  scoring  region  IR  , over  energy  bin  IE,  time  IT,  and  given  for 
all  IR,  IE,  IT,  for  each  statistical  aggregate  IA,  IA=1,  NOA,  where  NOA=NHIST/ 
NSTAT.  The  normalization  is  NSTAT  times  the  one  specified  by  the  SOUCAL  input. 
SAMCAR  calculates  the  sums 

NOA 

S ( IR, IE , IT)  = Z F ( IR, IE , IT , IA) 

IA=1 

NOA 

SQ ( IR , IE , IT)  = Z (F (IR, IE , IT, IA) ) 

IA=1 


* 

The  "effective"  volume  of  a point  detector  is  unity. 
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the  relative  statistical  errors 


cr(IR,  IE,  IT)  = 


NO  A 


S (IR,IE,IT) 


SQ (IR, IE, IT) 


NOA 


M 


1/2 


S (IR,  IE,  IT)\ 
NOA 


and  the  average  fluxes 


<MIR,IE,IT)  = 


S (IR, IE, IT) 


NHIST.V  (IR)*&:  (IE)  • AT  (IT) 


where  V(IR)  is  the  volume  of  scoring  region  IR,  £E( IE)  and  (S' (IT)  are  the 
widths  of  the  energy  and  time  bins. 

The  relative  statistical  errors  are  multiplied  by  100  and  are  quoted  in 
percent. 

6.3  Description  of  Input  Data  for  Program  SAMOUT 

Item  1 - Title  Card,  80  columns  of  Hollerith  (Format  20A4) 

Item  2-20  Available  Control  Switches  (Format  2012) 

KSWTCH(l)  = 0 Statistical  tape  NAGG  has  been  written 

= 1 statistical  tape  IAGG  must  be  processed  and 
NAGG  to  be  written 

KSWTCH(2)  = N N = number  of  response  functions  to  be  read. 

(Response  functions  of  unity  are  built  in  as 
a default  option; 

N=0  if  this  is  the  only  response  function. 

If  NAGG  has  been  made  on  a previous  run, 
this  switch  is  ignored.) 

KSWTCH (3)  = 0 Edit  all  problems. 

= 1 Omit  this  edit. 

KSWTCH (4)  = 0 Edit  all  problems  summed  over  time. 

= 1 Omit  this  edit. 


For  a point  detector,  the  volume  is  internally  set  to  unity. 
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KSWTCH (5)  = 0 Edit  all  problems  summed  with  response  functions. 

= 1 Omit  this  edit. 

KSWTCH(6)  = 0 Edit  all  problems  with  response  functions  summed 
over  time. 

= 1 Omit  this  edit. 

KSWTCH  ( 7 ) = N Number  of  problem  differences  required  for  this 
run. 

KSWTCH (8-11)  = 0 or  1 Same  as  for  KSWTCH (3-6)  except  here 
applied  to  required  problem  differences. 

KSWTCH (15-16)  = 0 No  intermediate  printout  (for  normal  usaqe) 

= 1 Print  intermediate  results 

If  KSWTCH (2)  = 0 or  KSWTCH (1)  = 0,  omit  next  item  3 and  go  to  item  4 


directly. 

Item  3 - If  KSWTCH (2)/0  and  KSWTCH (1)  = 1 only:  (Format  4E20.8) 
Read  prescribed  list  of  response  functions.  Each  new 
response  function  should  start  on  a new  card. 

If  KSWTCH (7)  = 0,  end  of  input;  otherwise  read  in  item  4. 

Item  4 - Prescribed  list  of  problem  differences  (Format  2014) 
e.g.  if  problem  differences  between  problem  2 and  4,  and 
between  problem  5 and  7 are  required,  this  card  will  read 
from  columns  1-8  the  following:  02040507. 

End  of  Data. 
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6.4  Tape  Utilization 

The  following  describes  the  function  of  each  tape  used  in  the  program. 
Tape  numbers  refer  to  Fortran  logical  numbers.  All  tapes  are  used  in  the 
binary  mode. 

Tape  1 

The  IAGG  tape,  the  statistical  aggregate  tape  from  the  previous 
program,  SAMCAR,  (denoted  as  TAPE  16  in  SAMCAR) . 

Tape  2 

A temporary  storage  tape  while  writing  the  NAGG  tape. 

Tape  3 

A temporary  storage  tape  while  writing  the  NAGG  tape. 

Tape  4 

The  NAGG  tape,  the  cumulative  statistical  tape  from  which  all 
editing  will  be  done. 
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7.  PROGRAM  SAMGAM:  SECONDARY  GAMMA  PRE-PROCESSOR 


7.1  General  Discussion 

Program  SAMGAM  is  the  data  processing  link  of  the  SAMCEP  system  that 
couples  a primary  neutron  transport  solution  to  its  corresponding  secondary 
gamma  transport  solution.  The  data  processing  involves  the  generation  or 
modification  of  three  data  files  required  for  the  secondary  gamma  transport: 

(1)  an  organized  (i.e.,  "BAND  FORMAT")  gamma  data  tape  (GODT) , generated  from 
a gamma  element  data  tape  (GEDT)  produced  by  SAM-X;  (2)  a perturbation  data 
tape  (PDT) , originally  produced  by  SAMIN , modified  by  SAMSAM,  and  again  modi- 
fied by  SAMGAM,  so  as  to  retain  only  concentration  perturbations  (the  only 
"neutron"  perturbations  relevant  to  the  secondary  gamma  transport) ; (3)  an 
external  source  tape  (EST)  produced  from  a neutron  interaction  file  (gener- 
ated in  the  primary  neutron  transport)  and  a gamma  production  data  tape  (pro- 
duced by  SAM-X) 

The  generation  of  the  GODT  and  PDT  is  performed  by  the  BANDG  and  PREP8 
routines,  respectively,  as  described  in  Section  7.3.  The  EST  is  generated  by 
the  SAMSOU  module  (i.e.,  set  of  routines  governed  by  SAMSOU) , which,  being 
the  principal  component  of  SAMGAM,  is  detailed  in  the  following  section.  Flow 
charts  for  the  SAMGAM  program  and  the  SAMSOU  module  are  given  in  Figures  7.1 
and  7 2,  respectively. 

7 2 The  SAMSOU  Module 

The  purpose  of  this  module  is  to  generate  from  a neutron  interaction  file 
(IF,  output  of  SAMCAR)  and  a gamma  production  data  tape  (GPDT) , the  external 
source  tape  (EST)  required  by  SAMCAR  for  secondary  gamma  transport. 

Each  neutron  interaction  represents  a non-elastic  event  for  each  of  up 
to  ten  correlated  problems.  The  information  distinguishing  any  one  problem 
from  all  the  others  is  a weight  which  is  the  cumulative  product  of  the  weights 
duo  to  the  relevant  perturbations  over  the  history  leading  to  the  event.  In 
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FIGURE  1.2  - Flow  Chart  for  SAMSOU  Module 


accordance  with  this,  a given  output  problem  may  also  be  made  the  subject 
of  perturbations  in  the  gamma  production  data  by  updating  the  secondary 
gamma  ray  weight  for  that  problem  for  each  relevant  perturbation. 

The  data  defining  a neutron  interaction  are  shown  in  Table  71.  Most 
of  the  items  are  self-explanatory,  but  item  number  10  may  be  clarified  by 
stating  that  KDLIV  contains  information  necessary  for  the  angle  reselection 
procedure  and  is  zero  if  there  are  no  detectors  (see  Appendix  F) 

As  depicted  in  Table  7.1,  the  data  defining  a secondary  gamma  ray  art  of 
the  same  type  and  are  recorded  on  the  source  tape  in  the  same  format  as  the 
interaction  data  for  a neutron  event.  The  cartesian  coordinates,  the  region 
number  and  the  time  of  the  source  gamma  are  identical  to  those  of  the  neutron 
event.  The  direction  and  energy  of  the  photon,  however,  are  selected  and  the 
various  problem  weights  must  be  updated  to  reflect  the  perturbation  data  and 
any  importance  sampling  used  in  selecting  the  photon  properties. 

The  general  flow  of  the  SAMSOU  module  is  shown  in  Figure  7.2  The  input 
consists  of  the  interaction  file,  the  gamma  production  data  tape,  the  organized 
gamma  ray  element  data  tape  (GODT) , and  card  input  specifying  the  perturbations 
t the  production  data  (if  any). 

ard  input  specifying  the  perturbations  and  the  detector  information  is 
id  by  tlie  subroutine  GIN  GIN  also  organizes  the  perturbation  and  energy  im- 
portance sampling  data  in  the  MISTER-SISTER  array  (not  to  be  confused  with 
identically  named  array  in  SAMCAR) . 

After  the  card  input  has  been  read,  the  next  task  of  SAMSOU  is  to  re- 
rganize  the  GPDT  so  that  only  those  elements  occurring  in  the  problem  are 
r presented  on  the  organized  gamma  production  data  tape  (OGPDT) . This  reorqani- 
at i in  i accomplished  by  a call  to  subroutine  SOUSEC,  which  utilizes  informa- 
• ion  from  the  gamma  ray  element  data  tape  (GODT)  to  select  the  elements  and 
■ ibli--.ii  the  order  of  occurrence. 
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The  next  task,  accomplished  through  a call  to  subroutine  ASORTT,  is  to 
sort  each  statistical  group  of  interactions  on  the  IF  so  that  the  element  order 
within  a group  corresponds  to  that  on  the  OGPDT.  The  output  of  this  operation 
is  a sorted  interaction  tape  (SIT) . 

The  main  calculational  routine  of  this  module  is  GAMGEN.  This  routine 
operates  on  a full  buffer  of  interactions  (35)  to  produce  the  required  number 
of  secondary  gamma  rays,  as  given  on  input.  A buffer  of  35  secondary  rays  is 
written  onto  the  source  tape  as  soon  as  it  is  filled. 

Two  types  of  perturbations  of  gamma  production  data  are  permitted: 

1)  perturbation  of  yield  (multiplicity)  data  for  specified  photons  from  speci- 
fied elements  interacting  with  neutrons  in  specified  energy  ranges  and  2)  pertur- 
bation of  angular  distribution  data  for  specified  photons  from  specified  elements 
in  collision  with  neutrons  in  specified  energy  ranges. 

In  addition  to  allowing  for  perturbation  of  the  gamma  production  data, 

SAMSOU  allows  for  the  selection  of  the  photon  from  an  alternative  yield  set, 
thus  providing  for  energy  importance  sampling.  The  alternative  yield  set  is 
input  as  a quadratic  function  of  gamma  energy  from  which  discrete  relative  yield 
values  are  obtained  for  the  given  photon  energies. 

The  secondary  gamma  ray  directions  which  are  recorded  on  the  source  tape 
are  subjected  to  the  reselection  procedure,  required  by  the  bounded  flux-at-a- 
point  (BFAP)  procedure,  in  GAMGEN.  Hence,  the  role  of  the  OGPDT  terminates  with 
the  completion  of  the  source  tape„ 

7,3  Description  of  Routines 

This  section  gives  a brief  description  of  the  routines  which  comprise  the 
AMGAM  program.  The  routines  are  presented  in  the  order  in  which  they  are  in- 
voked : 
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ROUTINE 


DESCRIPTION 


SAMGAM 


BANDG 


REED 


DSPLAY 


PREP8 


Sd  u .ouU 


This  is  the  main  program.  Its  principal  function  is 
to  direct  the  processing  of  the  three  data  files: 

GOD T,  PDT,  and  EST.  It  reads  one  input  card,  which 
specifies  whether  or  not  the  GODT  is  available  from 
a previous  execution  of  BANDG.  It  also  establishes 
the  length  of  blank  common  and  initializes  certain 
flags 

This  subroutine  is  called  by  SAMGAM  if  a GODT  has  not 
previously  been  prepared.  It  is  a simplified  version 
of  the  BAND  routine  of  program  SAMSAM.  The  GODT  is 
generated  in  one  energy  band 

This  utility  routine  provides  the  calling  routine 
(BANDG)  with  the  flexibility  of  effecting  an  unformatted 
read  of  a complete  array,  whose  dimension  is  an  argu- 
ment of  the  call.  This  obviates  the  use  of  an  implied 
loop  in  the  FORTRAN  supplied  READ  utility,  which  is 
slower. 

This  utility  routine  accepts  as  arguments  the  first 
member  of  a sequence  in  memory,  and  the  desired  number 
of  consecutive  entries  to  be  displayed.  It  then  prints 
out  the  specified  set  of  entries,  up  to  ten  per  line, 
with  an  internally  computed  cumulative  count  as  an 
eleventh  printed  integer.  It  automatically  discrimin- 
ates between  fixed  and  floating  point  entries,  and 
uses  either  112  or  E12.4  format,  respectively 

This  routine,  called  by  SAMGAM,  modifies  the  PDT  as 
follows:  every  entry  of  the  NTYP  array  (see  Section  3.1) 

is  set  equal  to  NTYP(l),  the  INBAND  array  (see  Section  4.  ) 
is  zeroed  out,  and  the  value  of  KEGEOM  (see  Section  4 
is  made  to  correspond  to  the  GODT.  The  PDT  is  then  over- 
written on  TAPE 8. 

The  governing  routine  for  the  external  source  generating 
module,  SAMSOU  is  called  from  SAMGAM  after  the  GODT  ha: 
been  produced  SAMSOU  reads  one  card  of  data  directly: 
NSTOP,  the  number  of  precursor  histories,  NSTAT  the  num- 
ber of  such  histories  per  statistical  group,  and  NGAM, 
the  desired  number  of  gamma  ray  histories.  The  rest  of 
the  card  input  is  then  read  by  GIN,  the  first  subroutine 
called  by  SAMSOU.  SAMSOU  then  calls  the  tape  organizing 
routines  SOUSEC  and  ASORTT.  After  ASORTT  has  returned 
the  total  number  of  interactions  on  the  SIT,  SAMSOU  cal- 
culates an  integer,  NFACT , and  a real  number,  SFACT  < 1.0, 
such  that  the  sum  represents  the  average  number  of  gammas 
to  be  produced  from  each  interaction.  Then  SAMSOU  reads 
in  buffers  of  35  interactions  from  the  SIT  and  calls 
GAMGEN  to  process  each  buffer 


ROUTINE 


DESCRIPTION 


GIN  This  subroutine,  the  first  called  from  SAMSOU, 

reads  the  rest  of  the  card  input,  and  locates 
data  for  the  perturbations  in  the  MISTER-SISTER 
array,  creating  pointers  as  it  reads  and  stores. 
The  data  read  by  GIN  also  includes  energy  im- 
portance sampling  parameters  and  detector  in- 
formation. The  detector  information  required 
is,  for  each  detector,  the  position  coordin- 
ates, the  radius  of  the  critical  sp’ ere,  the 
physical  region  and  the  overlap  index.  After 
all  of  the  data  have  been  read,  GIN  prints  it  out 
in  an  easily  read  format. 

SOUSEC  The  third  routine  called  by  SAMSOU,  SOUSEC  is 

responsible  for  organizing  the  gamma  production 
data  tape.  For  this,  S IUSEC  needs  the  header 
information  from  the  GODT,  specifically  the 
identification  of  elements  in  the  problem. 

SOUSEC  searches  the  GPDT  for  the  data  for  each 
element  and  writes  the  data  on  the  OGPDT  as 
two  logical  records  The  first  is  simply  the 
length  (LENGTH)  of  the  data  record.  The  organi- 
zation of  the  second,  or  data,  record  may  be 
found  in  Appendix  D of  the  SAM-CE  manual, 
starting  with  word  number  3 (the  first  two 
words  are  deleted  for  the  OGPDT) . The  pointers 
for  the  data  are  adjusted  to  reflect  the  known 
starting  location  in  the  MISTER  array,  which  is 
supplied  by  GIN.  SOUSEC  leaves  resident  ii 
memory  the  ordered  array  of  element  identifiers. 
The  only  other  output  is  the  OGPDT. 

ASORTT  The  function  of  this  routine  is  to  reorder  the 

interactions  within  each  statistical  group 
according  to  the  element  involved  in  the  i.tter 
action.  ASORTT  also  counts  the  number  of  inter 
actions  on  the  tape  (NINT)  for  use  in  computing 
the  number  of  gamma  rays  to  be  produced  for 
each  interaction.  The  input  to  ASORTT  is  the 
IF.  Other  required  information  is  the  ordered 
array  of  element  identifiers.  The  output  is 
the  sorted  interaction  tape  SIT. 


ROUTINE 


DESCRIPTION 


GAMGEN  GAMGEN  (IFLAG)  is  the  main  calculational  routine  of  this 

module.  This  routine  treats  a complete  buffer  of  35 
interactions,  producing  as  many  gamma  rays  per  interaction 
as  required  by  input.  For  each  gamma  ray  an  energy  (one 
of  a discrete  set)  and  a direction  must  be  selected.  The 
data  for  these  selections  are  read  into  the  MISTER-SISTER 
array  from  the  OGPDT  for  a given  element  and  remain  there 
until  the  interactions  for  that  element  are  exhausted. 

The  photon  energy  may  be  selected  from  the  true  yield 
data  or  from  an  alternative  set,  with  appropriate  weight- 
ing The  alternative  set  is  calculated  for  each  event 
from  a user-supplied  quadratic  function  of  photon  energy, 
the  same  for  all  neutron  energies.  Thus  some  energy  im- 
portance sampling  is  available. 

The  initial  direction  selection  is  based  on  the  data  from 
the  OGPDT  and  then,  if  detectors  are  present,  angle  re- 
selection with  the  indicated  re-weighting,  is  effected  bv 
a call  to  subroutine  ARG. 

After  the  energy  and  angle  selections  are  complete,  GAMGEN 
examines  each  perturbation  in  turn  to  determine  whether 
it  applies  to  this  element,  this  photon,  and  this  neutron 
energy.  The  additional  weight  factor  is  computed  for  each 
relevant  perturbation  and  multiplies  the  weight  of  each 
output  problem  to  which  the  perturbation  applies. 

GAMGEN  writes  buffers  of  35  gamma  rays  onto  the  source 
tape. 

Other  routines  called  by  GAMGEN  are  DIREC,  RANF,  SEEK,  and 
ARG,  which  calls  ARPREP.  These  are  described  elsewhere. 
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7.4  Card  Input  Formats 

This  section  describes  the  card  input  required  for  executing  the  SAMGAM 
program.  The  input  formats  are  presented  in  two  sections,  for  clarity,  but 
are  understood  to  represent  one  contiguous  data  record  The  two  sections  deal 
with  generation  of  the  GODT  and  the  EST,  in  that  order. 

7.4.1  Input  Formats  for  GODT  Generation 
Card  1;  Format  (215)  (read  by  SAMGAM) 

KEGEOM  Zero  (or  blank)  if  generation  of  GODT  ("banding')  is 

required;  otherwise,  the  value  of  KEGEOM,  as  printed 
out  in  the  previously  generated  GODT  run. 

I DBG  A positive  integer  value  will  cause  the  GODT  to  be 

displayed. 

If  KEGEOM  is  given  as  a positive  integer,  the  remainder  of  the  input  de- 
scribed in  this  sub-section  is  bypassed  (i.e.,  proceed  with  Section  7.4.2). 
Otherwise,  the  following  cards  are  read  by  BANDG. 

Card  2 : Format  (3110) 

NO  Arbitrary  run  identifier. 

NCOMP  Number  of  compositions. 

NG  Unity  for  gamma  cross  sections. 

For  each  of  the  NCOMP  compositions,  the  following  sets  of  cards  are  re^u  i. 
Card(s)  3:  Format  (110) 

NE  Number  of  elements  for  this  composition. 

For  each  of  the  NE  elements,  one  card  of  the  following  form  is  required: 
Card(s)  4:  Format  (2110,  El 5 6) 

IT  Arbitrary  counter. 

ID  Element  identifier  (ZZAAA) 

24  3 

CONt  Atomic  concentration  given  in  units  of  (10  atoms/rm  ) 

The  SAMGAM  input  continues  with  the  card  formats  described  in  the  section 


that  fol lows 


Input  Formats  for  Gamma  Ray  Source  Generation 
Card  1:  Format  (4110)  (read  by  SAMSOU) 

NSTOP  number  of  precursor  neutron  histories 

NSTAT  number  of  precursor  histories  per  statistical  group 

NGAM  total  number  of  gamma  histories  to  be  run;  default 

(if  left  blank)  leaves  decision  up  to  gamma  product- 

«* 

ion  data; 

NAGG  number  of  aggregates  to  be  used  in  a preliminary 

computation  of  the  multiplicative  factor  which  will 
produce  the  desired  NGAM  Default  (if  left  blank)  is  1. 
Card  2:  Format  (315) 

KK1  number  of  perturbations  of  type  1 

KK2  total  number  of  perturbations  (type  1 + type  2) 

Cards  3 through  8 are  needed  only  if  KK2>0. 

Card  3:  Format  (1415) 

Start  new  card  for  each  perturbation: 

IP  perturbation  type  (1  or  2) 

NPROB  number  of  output  problems  affected 

JPROB(l)  ordinal  number  of  first  affected  problem 

JPROB (NPROB)  ordinal  number  of  last  affected  problem 

The  rest  of  the  cards  are  grouped  by  perturbation  For  both  perturbation 

types : 

Card  4:  Format  (1415) 

ID  element  identifier 

NPH  number  of  photons  affected 

JPH(l)  ordinal  number  of  first  photon 

JPH(NPH)  ordinal  number  of  last  photon 
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Card  5: 


Format  (I5,2E15.5) 


NEP  number  of  neutron  energies 

EHIGH  highest  neutron  energy 

ELOW  lowest  neutron  energy 

Card  6:  Format  (5E15.5) 

(EN (I) ,1=1, NEP)  neutron  energies  for  which  perturbation  data 

are  tabulated.  EN (1) >EHIGH 
EN (NEP) < ELOW. 

For  perturbations  of  type  1: 

Card  7:  Format  (5E15.5)  (NEP  sets) 

(Y (1,J) ,J=1, NPH)  yields  for  EN(1)  for  each  photon 

(V (NEP, J) , J=1,NPH)  yields  for  EN (NEP)  for  each  photon 

NOTE:  Start  new  card  for  each  energy  point. 

'or  perturbations  of  type  2: 

Card  7:  Format  (1415) 

(NX (I) , 1=1, NEP)  number  of  CHI  values  for  each  energy  EN(I) 

Card  8:  Format  (5E15.5)  (NEP  sets) 

iCHI (1, J) , J=1,NX (1) ) CHI  values  for  EN(1) 

I 

(CHI (NEP,J) ,J=1,NX(NEP)  CHI  values  for  EN (NEP) 

NOTH:  Endpoint  values  (0.0  and  1.0)  are  not  to  be  given 

in  the  CHI-sets. 
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Card  9:  Format  (3E15.5) 

A 1 Coefficients  in  the  quadratic  formula 

B *•  for  alternative  yield  data:  (Y(EI)=AE^  + BE^.  + C. 

C Can  be  used  to  obtain  some  energy  weighting.  For 

no  energy  importance  enter  (0.0, 0.0, 1.0) . 

Card  10:  Format  (15) 

NDET  Number  of  detectors. 

For  each  of  the  NDET  detectors: 

Card  11:  Format  (4E15.5,2I5) 

XAD  | 

YAD  x,  y and  z coordinates  for  this  detector. 

\ 

ZAD  ' 

CRAD  radius  of  the  critical  sphere  for  this  detector,  output 

of  the  neutron  problem 

IRDET  physical  region  of  this  detector,  output  of  the  neutron 
problem. 

I.AP  signal  indicating  either  the  absence  (0)  or  the  presence 

(1)  of  overlays  of  the  critical  sphere  of  this  detector 
with  that  of  any  other  detector.  Also  an  output  of  the 
neutron  problem. 
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Table  7.1  The  Twenty-Four  Word  Array  Defining  a Non-elastic  Event 
and  the  Equivalent  Array  Defining  a Secondary  Source 

Description 


Word (s) 

Name 

Non-elastic  Event 

Secondary  Gar. 

1-3 

XB 

Cartesian  Coordin- 
ates of  event 

= 

Cartesian 
ates  of  Gal 

4-6 

WB 

Direction  of  pre- 
cursor primary 
particle 

Direction  s« 
for  Gamma 

7 

E 

Energy  of  precursor 
primary  particle 

Energy  oi 
Gamma 

8 

IR 

Region  number 

- 

Reuioii  num  • 

9 

T 

Time  of  event 

= 

Time  of  ori 

10 

IDET 

200000*KDLIV 
+SIGN (KDLIV) *IATWT 

= 

Same 

11 

F 

Ignored 

Ignored 

12 

NHIST 

Primary  Particle 
History  number 

Same 

13 

W 

Ignored 

Ignored 

14 

J12345 

10 

10 

15 

WT  (1) 

Weight  for  Problem 

1 

Updated  W<  ■ : 
for  Probler 

16 

WT  (2) 

Weight  for  Problem 

2 

Updated  Vo  i 
for  Pi  obi  i.. 

24 

WT (10) 

Weight  for  Problem 

10 

Updated  W<  i i!v 

Problem  l'1 
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7 5 Organization  of 

the  MISTER-SISTER 

Array 

(This  array  is 

independent  of  the 

identically  named  array  in  the  SAMCAR 

program) . 

LOCATION 

CONTENTS 

REMARKS 

1 

LN1 

Location  of  NPROB  for  perturbation  1 

KK 

LNKK 

Location  of  NPROB  for  perturbation  KK 

KK2 

LNKK2 

Location  of  NPROB  for  perturbation  KK2 

KK2+1 

LOCO 

Location  of  last  output  problem  number' 
for  last  perturbation 

KK2+2=LNl 

KK2+3 

NPROB 

JPROB(l) 

(NPROB  for  first  perturbation) 
first  output  problem  number  for 
first  perturbation 

KK2+2+I 

JPROB(I) 

Ith  output  problem  number  for 
first  perturbation 

KK2+2+NPROB 

JPROB (NPROB) 

last  output  problem  number  for 
first  perturbation 

(Bracketed  sequence  repeated  for  each  perturbation) 


• 

• 

LOCD+1 

LOCI 

location  of  first  data 

word 

for 

. 

. 

first  perturbation 

LOCCH-I 

LOCI 

location  of  first  data 

word 

for 

• 

Ith  perturbation 

LOCD+KK2 

LOCKK2 

location  of  first  data 
last  perturbation 

word 

for 

1.14 


L0CD+KK2+1 


LOCIMP 


location  of  first  coefficient  of 
importance  function 


(Following  sequence  for  each  perturbation  in  order) 


LOCATION 

CONTENTS 

REMARKS 

LOCI+1 

ID 

element  ID  for  Ith  perturbation 

I.OCI+2 

NPH 

number  of  photons  affected 

LOCI+3 

JPH(l) 

ordinal  number-first  photon 

L0CI+2+J 

JPH(J) 

• 

ordinal  number- Jth  photc  . 

L0CI+2+NPH 

• 

JPH(NPH) 

ordinal  number-last  phot  n 

LOCI+NPH+3 

NEP 

number  of  energies  in  neuti  >i 

LOCI+NPH+4 

EHIGH 

high  energy  limit 

LOCI+NPH+5 

ELOW 

low  energy  limit 

LOCI+NPH+6 

— 

blank 

LOCI+NPH+7 

EN  (1) 

first  energy,  neutron  table 

I CI+NPH+6+J 

EN  ( J ) 

Jth  energy,  neutron  table 

• 

LOC I +NPH+6+NEP 

EN (NEP) 

• 

last  energy,  neutron  table 

NOTE:  the  sequences  for  the  two  types  of  perturbation  diverge  here 

i • t LOCIP-LOCI +NPH+6+NEP . 

For  Type  1 : 

LOOIP+1  Y (1 , 1)  yield  for  photon  JPH(l)  at  £N  i i ' 

IjOCIP+NPH  Y (1  ,NPH)  yield  for  photon  JPH(NPH)  at  EN  1 1 , 

• • • 

0CIP+N£P*NPH+1  Y (NEP , 1 ) yield  for  photon  JPH(l)  at  t !.  ■ 

yield  for  ptioton  JPH(NTH)  it  I 

1« 


tl’t  ( NEP-t  1 ) *NPH 


Y (NEP.NPH) 


For  type  2: 


LOCATION 

CONTENTS 

REMARKS 

LOCIP+1 

NX  ( 1 ) 

number  of  CHI  values  at  £N(1) 

LOCIP+NEP 

NX (NEP) 

number  of  CHI  values  at  EN(NEP) 

LOCIP+NEP+1 

CHI (1,1) 

CHI-1  at  energy  EN(10) 

LOCIP+NEP+NX ( 1 ) 

CHI ( 1 ,NX ( 1 ) ) 

CHI-NX(l)  at  energy  EN(1) 

• 

• 

• N 

• 

• 

\ 

NEP— 1 

LOCIP+NEP+  E NX  (J) +1 

CHI (NEP , 1 ) 

CHI-l  at  energy  EN(NEP) 

J=1 

* 

• 

* 

O 

• 

NEP 

LOCIP+NEP+  l NX ( J) 

CHI (NEP, NX (NEP) ) 

CHI-NX (NEP)  at  EN(NEP) 

J=1 

LOCIMP 

A 

coefficient  of  square  of  photon 
energy 

LOCIMP+1 

B 

coefficient  f photon  ercrgy 

LOCIMP+2  =LOCEL 

C 

constant 

- .:A,B,C  '.re  used  to 

compute  alternate 

yield  data  (relative  yields)  . 

LOCEL+1 

ID 

element  identifier  - start  of 
gamma  production  data 

7.6  Tape  Utilization 

The  following  describes  the  function  of  each  tape  used  in  the  proat  : 
Tape  numbers  refer  to  FORTRAN  logical  numbers.  All  tapes  are  used  in  t!  • 
binary  mode. 

Tape  8 

This  is  the  perturbation  data  tape  (PDT) , generated  by  SAMIN,  then  mo- 
fied  by  SAMSAM  and  by  SAMGAM  (retaining  only  concentration  perturbatio 
if  any) . 

Tape  10 

This  is  the  organized  data  tape  for  gammas  (GODT)  generated  by  SAMo  • 
Tape  11 

This  is  the  gamma  element  data  tape  (GEDT)  generated  by  SAM-X,  from 
which  the  GODT  is  generated. 

Tape  12 

This  is  the  gamma  production  data  tape  (GPDT)  generated  by  SAM-X. 

Tape  1 4 

This  is  the  neutron  interaction  file  (IF)  generated  in  the  primary 
neutron  run  by  SAMCAR 
Tape  15 

This  is  the  external  source  tape  (EST)  produced  by  SAMGAM. 
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APPENDIX  A 


THE  COMBINATORIAL  GEOMETRY  TECHNIQUE  FOR  THE  DESCRIPTION 
AND  COMPUTER  PROCESSING  OF  COMPLEX  THREE-DIMENSIONAL  OBJECTS 

A I Introduction 

Combinatorial  Geometry  is  essentially  a technique  for  representing,  in 
a computer,  a mathematical  model  of  a three-dimensional  geometric  configura- 
tion. Once  in  the  computer,  the  configuration  can  be  analyzed  in  mam  di  : 
ent  ways  by  ray-tracing  techniques.  For  example,  quantities  such  as  vi  , 
surface  areas,  object  boundaries,  line  of  sight  distances,  etc.  are  readi 
determined.  Regardless  of  the  application,  however,  the  basic  concepts  em- 
ployed are  the  same.  A discussion  of  these  concepts  can  logically  be  broker 
down  into  two  topics.  That  is,  geometry  description  and  ray-tracing,  which 
are  discussed  separately  below. 

A. 2 Description  of  the  Combinatorial  Geometry  Technique 

In  order  to  perform  computer  studies  concerning  a complex  three-dimen- 
sional object  one  must  first  be  able  to  prepare  a mathematical  model  of  the 
object  and  its  environment.  The  combinatorial  geometry  technique  has  bet  r. 
developed  to  permit  a model  to  be  produced  which  is  both  accurate  and  suital . 
for  a ray-tracing  analysis  program.  The  latter  feature  is  important  since 
nuclear  radiation  analysis  by  Monte  Carlo  involves  the  tracing  of  rays  through 
geometrical  models 

In  effect  the  geometric  description  subdivides  the  problem  space  into 
unique  regions.  This  is  achieved  through  the  use  of  ten  specific  geometric 
bodies  (closed  surfaces)  and  the  orderly  identification  of  the  combination  of 
these  bodies  to  define  a region  (space  volume).  The  bodies  which  will  he 
!;■. cussed  further  in  INPUT  Prex>aration  (Section  A. 2 2)  are  as  follows: 


1.  Rectangular  Parallelepiped  (RPP) 

2.  Box 

3.  Sphere 

4.  Right  Circular  Cylinder 

5.  Right  Elliptical  Cylinder 

6.  Truncated  Right  Angle  Cone 

7.  Ellipsoid  of  Revolution 

8.  Right  Angle  Wedge 

9.  Arbitrary  Convex  Polyhedron  of  four,  five  or  six 
sides  (each  side  having  three  or  four  vertices) 

10.  Truncated  Elliptical  Cone 

Except  for  the  RPP's,  all  bodies  may  be  arbitrarily  oriented  with  respect 
to  the  x,  y,  z coordinate  axes  used  to  determine  the  space.  It  should  be 
noted  that  the  sides  of  an  RPP  must  be  parallel  to  the  coordinate  axes. 

A. 2,1  Region  Description  Technique 

The  basic  technique  for  the  description  of  the  geometry  consists  rf  de- 
fining the  location  and  shape  of  the  various  physical  regions  (wall,  equ i 
ment,  etc.)  in  terms  of  the  intersections  and  unions  of  the  volnr.es  co- 
in a set  of  simple  bodies.  A special  operator  notation  involving  the  syi 
l+j,  (-)  , and  (OR)  is  used  to  describe  the  intersections  and  unions.  Thes 
symbols  are  used  by  the  program  to  construct  tables  used  in  the  ray-tracing 
portion  of  the  problem. 

If  a body  appears  in  a region  description  with  a (+)  operator,  it  means 
that  the  region  beinq  described  is  wholly  contained  in  the  body. 

If  a body  appears  in  a region  description  with  a (-)  operator,  it  means 
that  the  region  being  described  is  wholly  outside  the  body. 
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The  (OR)  operator  is  used  to  form  regions  as  unions  of  subregions, 
where  each  subregion  is  defined  in  terms  of  one  or  more  bodies,  using  (+) 
or  (-)  as  described  above.  Then  a point  is  in  the  region  if  it  is  in  any 
subregion. 

The  technique  of  describing  a physical  region  is  best  illustrated  by 
an  example.  Consider  an  object  composed  of  a sphere  into  which  is  inserted 
a cylinder.  This  is  shown  in  cross  section  in  Figure  A. 2. 1(a). 

To  describe  the  object,  we  take  a spherical  body  penetrated  by  a cl. 
drical  body  (Figure  A02.1(b)).  Each  body  is  numbered.  Consider  the  spht 
as  body  No.  1 and  cylinder  as  body  No.  2.  If  the  materials  in  the  sphe; 
and  cylinder  are  the  same,  then  they  can  be  considered  as  one  physical  region, 
say  region  100  (Figure  A.2.1(c)). 

The  description  of  region  100  would  be 
100  = (OR  1)  (OR  2). 

This  means  that  a point  is  in  region  100  if  it  is  either  inside  body  1 or 
inside  body  2. 

If  different  materials  are  used  in  the  sphere  and  cylinder,  then  the 
sphere  with  a cylindrical  hole  in  it  would  be  given  a different  reqion  numb, 
(say  200)  from  that  of  the  cylinder  (300). 


I'll 


The  description  of  region  200  would  be  (Figure  A. 2. 1(d)) 


200  = (+1)  (-2). 

This  means  that  points  in  region  200  are  all  those  points  inside  body  1 which 
are  not  inside  body  2. 

The  description  of  region  300  is  simple  (Figure  A.2.1(e)) 

300  = (+2) 

That  is,  all  points  in  region  300  lie  inside  body  2. 

This  technique,  of  course,  can  be  applied  to  combinations  of  more  th, 
two  bodies  and  such  region  descriptions  could  conceivably  contain  a long 
of  (+) , (-)  and  (OR)  operators.  The  important  thing  to  remember  is  that 

spatial  point  in  the  geometry  must  be  located  in  one  and  only  one  region, 
examples  are  given  in  Section  A. 2. 2. 2. 

The  user  of  the  program  will  specify  the  geometry  by  establishing  two  tal 
The  first  table  will  describe  the  type  and  location  of  the  set  of  bodies  use  : 
in  the  geometrical  description.  The  second  table  will  identify  the  physical 
region  in  terms  of  these  bodies.  The  computer  program  processes  these  tables 
to  put  the  data  in  the  form  required  for  ray  tracing.  All  of  the  space  must  t • 
divided  into  regions,  and  once  again,  no  point  may  be  in  more  than  one  regjw. 

A. . '.2  Input  Preparation 

A. . . . 1 Description  of  Input  Parameters 

The  information  required  to  specify  each  type  of  body  is  as  follows. 
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These  bodies  are  used  for  gross  subdivisions  of  the 
geometry  and  must  have  bounding  surfaces  parallel  to  the 
coordinate  axes.  Specify  the  maximum  and  minimum  values 
of  the  x,  y,  and  z coordinates  which  bound  the  parallelepiped. 


Box  (BOX) 


Specify  the  vertex  V at  one  of  the  corners  by  giving  its 
(x,y,z)  coordinates.  Specify  a set  of  three  mutually 
perpendicular  vectors,  H^,  representing  the  height,  width 
and  length  of  the  box,  respectively.  That  is,  the  x,y, 
and  z components  of  the  height,  width,  and  length  vectors 
are  given. 


Sphere  (SPH) 

Specify  the  vertex  V at  the  center  and  the  scalar,  R, 


denoting  the  radius. 


Truncated  Right  Angle  Cone  (TRC) 


Specify  a vertex  V at  the  center  of  the  lower  base,  the 
height  vector,  H^,  expressed  in  terms  of  its  x,  y,  2 com- 
ponents, and  two  scalars,  R^  and  R^,  denoting  the  radii 
of  the  lower  and  upper  bases. 


a. 


Ellipsoid  (ELL) 

Specify  two  vertices,  V ^ , denoting  the  coordinates  of  the 
foci  and  scalar,  R,  denoting  the  length  of  the  major  axis. 
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Right  Angle  Wedge  (RAW) 

Same  input  as  for  the  boxes.  However,  and  describe 
the  two  legs  of  the  right  triangle  of  the  wedge. 


Arbitrary  Polyhedron  (ARB) 

Assign  an  ordinal  number  (1  to  8)  to  each  vertex.  For  each 
vertex,  give  the  x,y,z  coordinates.  For  each  side  of  the  figure 
list  the  ordinal  vertex  numbers.  The  vertices  and  side  descrip- 


tions may  be  given  in  any  order.  An  example  is  given  later. 


Truncated  Elliptical  Cone  (TEC) 

Specify  the  coordinates  of  the  vertex  V at  the  center 
of  the  larger  ellipse;  the  x,y,  and  z components  of 
height  vector  II;  the  components  of  normal  vector,  N, 
directed  inward  at  V;  the  components  of  direction  vector, 
A,  along  major  axis;  the  semi-major  and  semi-minor  axes 
of  larger  ellipse,  and  R^,  respectively;  the  ratio,  P, 
of  the  larger  to  the  smaller  ellipse  axis.  Note  that 
direction  vectors  N and  A are  normalized  internally  (aftei 


input  printout) . 


A, 2. 2. 2 Examples  of  Region  Descriptions 


Some  representative  geometries  and  their  input  descriptions  are  shown 
below. 

Example  1 - Two  Spheres  Within  an  RPP  (See  Figure  2.2) 

The  body  input  table  is  shown  below. 

TABLE  I - BODY  INPUT  DESCRIPTION 


Body  Type  of  Data  Required 

1 List  the  six  bounding  coordinate  values 


(X.,X  ,Y.,Y  ,Z.fZ  ) 

min  max  nun  max  min  max 

2 List  the  vertex  and  radius  of  sphere  2 

List  the  vertex  and  radius  of  sphere  3 
One  possible  region  input  table  is  shown  below. 

TABLE  II  - REGION  DESCRIPTION 


Reg  ion 
100 


200 


300 


400 


500 


I nput 

(+1)  (-2)  (-3)  (Region  100  is  composed  of  all  points 

interior  to  RPP  No.  1 and  exterior  to  spheres  2 and  3) 

(+3)  (-2)  (Region  200  is  composed  of  all  points  interior 

to  sphere  3 and  exterior  to  sphere  2) 

(+2)  (+3)  (Region  300  is  composed  of  all  points  which  are 
in  sphere  2 and  are  also  in  sphere  3) 

(+2)  (-3)  (Region  400  is  composed  of  all  points  interior 
to  sphere  2 and  exterior  to  sphere  3) 

(OR  2)  (OR  3)  (If  desired,  one  region,  the  total  of 
regions  200,  300,  and  400,  can  be  defined  as  region  SQ0) . 


Example  2 - Cylinder  Divided  into  Two  Regions  by  a Box  and  with  a 


Sphere  at  One  End  (See  Figure  A. 2. 3) 

TABLE  I - BODY  INPUT  DESCRIPTION 


1 


5 


3 


-1 


The  req ion 

Region 

100 


200 


300 


400 


Type  of  Data  Required 

List  the  six  bounding  coordinates  of  the  RPP 
List  the  vertex,  radius,  and  height  vector  of 
cylinder 

List  center  and  radius  of  sphere 
List  coordinates  of  one  corner  and  components 
of  three  vectors  representing  sides  of  box. 
input  is  as  follows. 

TABLE  II  - REGION  DESCRIPTION 
Input 

(+1)  (-2)  (-3)  (All  points  interior  to  the  RPP  and 

exterior  to  the  cylinder  and  sphere.  Note  that 
region  100  includes  all  of  the  space  contained  inside 
body  4,  except  that  portion  inside  cylinder  2.  This 
space  can  be  assigned  a special  region  number,  if 
desired.  If,  as  in  this  example,  it  is  not  desired, 
it  is  not  necessary) . 

(+2)  (-4)  (All  points  interior  to  the  cylinder,  and 
outside  the  box) . 

(+3)  (-2)  (All  points  interior  to  the  sphere  and 

external  to  the  cylinder) . 

(+2)  (+4)  (All  points  interior  to  the  cylinder  and 

also  inside  the  box) . 
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Example  3 - Multiple  Region  Capability  - Cylinder  Containin' 


Two  Spheres,  ftll  Inside  an  RPP  (See  Figure  A. 2.4) 

TABLE  I - BODY  INPUT  DESCRIPTION 
Type  of  Data  Required 
List  RPP  data 
List  cylinder  input 
List  sphere  input 
List  sphere  input 
TABLE  II  - REGION  DESCRIPTION 
Input 
(+1)  (-2) 

(OR  3)  (OR  4)  (All  points  interior  to  3 or  4) 
(+2)  (-3)  (-4)  (All  points  in  the  cylinder  but 

not  in  the  spheres) . 


Body 

1 

2 

3 

4 

Region 

100 

200 

300 


1 *34 


A 2.2  3 Card  Input  Formats 

The  following  punched  cards  are  needed  to  describe  the  geometry  and  must 
appear  in  the  order  in  which  they  are  described  below. 

1.  Body  Cards 

The  computer  assigns  to  each  body  an  ordinal  number  which  depends  on  the 
order  in  which  the  body  cards  are  read  in.  Therefore,  it  is  most  important  that 
the  card  sequence  match  the  numbering  sequence  used  in  the  region  descriptions. 
Note  hat  no  gaps  may  be  left  in  the  body  numbering  sequence. 

Ten  different  body  types  may  be  employed.  The  standard  format  for  each 
is  follows: 


Plank 

Three-letter  body  identifier 
Blank 

Four  characters  or  arbitrary  integer  data 
Divided  into  six  floating  point  fields  of  10 
columns  each.  Body  dimensions  t giv  i. 

’.1  describes  the  input  required  for  * * h bo;  . T qua  i .<■  \ 

■ c'  lined  in  Section  A.  2. 2. 

■ tha  the  last  card  of  the  body  data  must  be  EM'  punc.ied  m coiu,..n: 


':'h ; ■;  is  r.he  signal  that  all  body  data  have  been  treated. 


INPUT  REQUIRED  FOP.  EACH  BODY  TYPE 


4-4 

o 


'O 

0) 

T3 

<D 

V 

z 


ui 

G 

U 


CM  CM 


4-4  4-4 

O 0 

rH  CM  .H 


CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

in 

in 

in 

in 

CO 

ro 

CO 

4-4 

4M 

4-4 

4-4 

4-4 

4-4 

4-4 

4-4 

4-4 

4-4 

4-4 

4-4 

44 

4-4 

4-4 

4-4 

4-4 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

rH 

CM 

i — 1 

CM 

iH 

CM 

rH 

CM 

rH 

CM 

rH 

CM 

ro 

<3* 

«—4 

CM 

CO 

CM  CM 
rH  CO 

X X 


CM 

X 


N 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

CM 

rH 

CO 

CM 

X 

X 

> l 

X f 

X 

X 

> 

CM  CM  N 
^ CO 
> > > 


CM  CM 
X < 


O 

I 


X 


>1 
> CM 

X X 


> I 


>1 

X 


.h  co 
X X 


^ ^ ^ 
cm  sr  vD  co 
> > > > 


0> 

CP 

CO 

x 

cp 

c 

•H 

8 


X < 


</} 

in 

X 

G 

i 

03 

X 

X 

iH 

£ 

rH 

co 

3 

*3* 

>« 

X 

X 

rH 

0 

u 

TJ 

O 

G 

C 

<TJ 

i 

•H 

CM 

o 

iH 

£ 

CM 

r j 

CO 

> 

> 

X 

0 
cn 

1 


0 

cm 

1 


>1 
>,  CM 
> X 


X 

X CM 
> X 


CM 

> 


> 


X 

> 


X 

X 


N 

> 


>1 

> 


X 

> x 


X 

X CM 
X X 


>1 
>1  rH 

> x 


X 

X rH 

> X 


X 

CM 

> 


N N 

N H iH 

> x 


> I 


> I 


> X 


X 

X 


CM 

> 


X X 


>1  CM 

> x 


X ' 
> X 


>1 
>,  CM 

> X 


X 

X CM 
> X 


X X X X 
CM  vD  CO 
> > > > 


N CM  CM  CM  CM 

N CM  H fO  ID  Is 

> X > > > > 


>*  >1  >i  >* 

h ci  in  p 

> > > > 


X X X X 
h ci  in  p 
> > > > 


0 

4-4 

a; 

P 

O 

z 

CD 

a> 

to 


w 

c 

c 

■H 

P 

X 

•H 

u 

u 

1/5 

<D 

Q 

CD 

U 

CO 

X 


X X 
X < 


CM  N 

> Z X 


>,  >,  CV 

> z x 


xx  - 
> z X 


X 

t- 


G 

Or 

X 

<D 

X 

o 

X 

4J 

4-> 

2 

X 

CO 

y 


U 

2 


u 


£ 


o 


TJ 

$ 

G 

U 

a> 

X 

rfl 

•H 

c 

H 

rH 

4J 

4-1 

0 

G 

X 

3 

X 

O G 

u 

CD 

u 

*1—4 

O 

rH 

4 

G 

G 

rH 

G 

TJ  -H 

n 

CP 

•H 

a> 

rH 

<D 

•H  G 

a> 

G 

u 

T> 

U 

TJ 

0 3 

G 

< 

03 

c 

c 

C/1  rH 

40 

a> 

G 

4J 

•H 

•G 

•H 

X 0 

u 

P 

CP 

G 

03 

-C 

rH 

.c 

1 

•H  > 

c 

rC 

TJ 

c?  X 

-C 

CP 

>1 

CP 

rH  <D 

3 

CP 

C 

X 

•H 

u 

■H 

U 

rH  X 

G 

•H 

2 

X 

X 

CO 

X 

X 

W 

E-* 

X 

c 

o 

>,  g 

G "0 


p 

X 


t; 

c 

* 

Ctf 

u 

c 

p 


lri7 


Note:  Each  of  the  six  faces  of  an  ARB  are  described  by  a four-digit  number 

giving  the  number  of  the  four  vertex  points  at  the  corners.  The  order  of 


specification  of  the  i>,ar  points  is  completely  arbitrary.  The  point  specifi- 
cation format  is  6(E10.3)  starting  in  Column  11.  An  example  is  shown  below. 

7 

y 

face  1 2 3 4 5 6 

PTS  , 1653.  3548.  4278.  1762.  1243.  5678. 

: 

Figure  A. 2.  - Example  of  Arbitrary  Polyhedron 

ii  the  number  of  faces  is  less  than  six,  the  remaini  iq  fa. 

ro,  and  must  appear  at  the  end  of  the  list. 

face  as  three  vertices,  the  omitted  position  may  i . i h.  . 

i the  other  vertices. 

t vortices  must  always  be  supplied.  Those  that  do  rot  appear  in 
Jr  are  ignored. 
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2.  Region  Cards 


Each  region  must  be  numbered  and  described  by  a logical  combination  of 
the  bodies  which  make  up  that  region.  Use  as  many  cards  as  necessary  to  de- 
scribe each  region  and  begin  each  region  on  a new  card.  The  input  format, 
described  below,  is  (2x,  A3,  5x,9  (A2,  15)). 


Columns 


1-2 

3-5 

6-10 

11-73 


Input 
Blank 

Arbitrary  hollerith  data 
Blank 

Divided  into  nine  fields,  of  7 columns  each.  The 
first  two  columns  of  each  field  are  reserved  for 
the  OR  operator.  The  third  column  is  for  the  (+) 
or  (-)  operator.  The  last  four  columns  are  for  the 
body  number. 

Use  as  many  cards  of  the  above  type  as  needed  to  complete  a region 
description,  but  leave  Columns  1-10  blank  on  all  continuation  cards 

The  last  region  description  card  must  be  followed  by  a card  contains  F 
in  Columns  3,  4,  and  5.  This  informs  the  code  that  all  regions  have  been  ■.  i 


scribed 
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APPENDIX  B 


ORGANIZATION  OF  THE  MASTER- ASTER  ARRAY 

The  first  part  of  the  MASTER-array  contains  the  unperturbed  cross  sec- 
tion data  in  the  same  format  as  the  output  of  SAM-X  (processor  of  ENDF/B  cross 
section  data) , the  details  of  which  are  referred  to  in  a separate  report  of 
processor  SAM-X1.  The  remainder  of  the  MASTER-array  is  described  below. 

The  capitalized  name  in  each  box  refers  to  a "pointer".  These  "pointti 
are  used  to  locate  sections  of  data  with  a minimum  of  calculation. 


Comments 


The  scores  are  accumulated  and  stored 
by  CARLO  and  FLUP.  The  scores  are 
stored  as  a function  of  energy  and 
scoring  region. 

This  section  of  the  array  contain: 
region-dependent  parameters  stored  two 
computer  words  per  region.  The  actual 
data  are  discussed  in  the  UNFF  se<  t;  01 
of  this  report.  The  data  art  r<  n 
and  stored  by  INPUTD. 


let 


PBSCKDIiG  tfAGfc.  bLANK-HOT  iTILHiD 


Contents 


Comments 


LBIRTH 

(4) 

The  number  of 

j 

births  per  region 

LDEATH 

(5) 

The  number  of 

j 

deaths  per  region 

j 

LESCAP 

(6) 

The  number  of 

escapes  per 

region 

LNDEG 

The  number  of 
degrades  per 
region 


LNABS 

*8)  The  number  of 

absorptions  per 
reg ion 


>r  it-am  •>  (3)  through  (8)  the  starting  locations  are  compi  ‘ e ; n IN1."' 
i”  counts  ar;  accumulated  in  CARLO  and  are  printed  by  TALLY. 


! 

LRA  V 

j 

(9) 

Region 

weights 

LREW 

| 

(Id) 

Region 

energy 

weight 

i 

sets 

' 

The  actual  region  weights  to  be  us-  1 
for  region  importance.  The  weights 
are  read  in  by  INPUTD. 

The  energy  weight  sets  for  imports- 
sampling.  This  section  exists  if 
energy  importance  is  used  in  the  pic' 
lem.  The  data  are  read  in  by  iNPtrrr . 


U>2 


Contents 


Comments 


The  aiming  angles  for  angular  importance. 

Three  words  per  angle  denoting  direction 
cosines  are  stored.  The  data  are  read  by 
INPUTD.  The  array  exists  only  if  angular 
importance  is  used. 

The  angular  weight  sets  for  angular  impor- 
tance. The  array  exists  only  if  angular 
importance  is  used.  The  data  are  read  by 
INPUTD. 

The  energy,  position,  and  direction  data  : 
the  source  distribution.  The  data  ar<  rea 
in  by  SOUCAL. 

This  section  uses  up  all  available  rooi 
the  MASTER-array . Supergroup  latents  ar 
stored  here.  Tapes  will  be  used  for  latent 
storage  if  the  available  room  is  insufficient. 
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APPENDIX  C 

ORGANIZATION  OF  THE  MISTER-SISTER  ARRAY 

The  following  sections  describe  the  organization  and  contents  of  the 
MISTER-SISTER  array  in  which  all  the  perturbation  data  and  necessary  inform- 
ation ab<5ut  them  reside. 

Cl  - Storage  Description  for  Problem  Description 


Contents 


Comments 

LOC's  are  location  pointers  which 
give  the  location  which  contains 
the  number  of  problems  (NPROB) 
affected  by  a specific  perturbafi'  i. 
(IP)  . 


NPER  = NTYP(IO)  = total  number  of 
perturbations 


LOCLIP  = last  address  of  problem 
description 


NPROB  words 


16S 


f riiiiCKDXi G cLAllK^iJOT  rlL-^D 

I m 


(All  notations  refer  to  text  of  Section  3.2) 


Contents  Comments 


LOCLIP+1 


LOCLIP+ 
n ER  + 1 

is  the  actual  perturbation  Hat  i 
classified  by  perturbation  t’l 


LOCIP 
for  IP=1 


LOCIP 
for  IP=2 


LOCIP 

for  IP=NPER 


BLANK 


LOCIP's  are  location  pointer 
which  give  the  location  where 
perturbation  data  of  a specif i 
perturbation  (IP)  begin. 


Following  location  (LOCLIP+NPEK+ I t 


TYPE  1 - Composition  (Present  if  NTYP ( i )>  0) 
Contents  Corrai . nts 


these  repeat  NTYP(1)-1  times 
if  NTYP(1)>1. 


TYPE  2 - Whole  Element  Pert  .rbation 


Contents 


Comments 


LOCIP+O 
+ 1 

+ 2 
+ 3 
+4 


+ 5 


LOC (LOCELM) 


LOC (LOCELP) 


EH 


EL 


LOSAM  for 
CHI -table 


LOSAM 

for  ENN-table 


Other  perturba- 
tions of  TYPE  2 
repeat  the  above 
layout 


LOC (LOCELM) =location  of  LOCELM 
in  MASTER-array 

LOCELM= location  where  the  un;  < r- 
turbed  data  of  an  element  begin 
in  MASTER-array. 


LOC (LOCELP ) -location  of  LO  RLi 
in  MASTER-array. 


LOCELP=location  where  the  pertu: 
element  data  begin  in  MASTER-array. 
(Only  exists  for  perturbation  type  . - 

LOSAM  = location  of  sampling  table. 


TYPE  3,  4 and  5 - Total,  Scattering  or  Inelastic  Scattering 
Cross  Section  Perturbation 


Contents 

LOCIP+O 

LOC (LOCEJM) 

+ 1 

NEP 

+ 2 

EH 

+ 3 

EL 

+4 

FMULT 

+ 5 

LOOPS 

E.  >EH 

Comments 


EMULT  is  significant  only  if 
NEP<2;  LOCIP+5  to  L0CIP+5+NEP 
exist  only  if  NEP^-2 


TYPE  3,  4 and  5 - continued 


Contents 


Comments 


a (Ex) 


a (E2) 


0(IW 


Other  perturbations 
of  TYPE  3 repeat  the 
above  layout 


Data  for  perturba- 
tions of  TYPE  4 


Data  for  perturba- 
tions of  TYPE  5 


Cross  section  data  at  given  energies 


TYPE  6 - Angular  Distribution  of  Elastic  Scattering 

Contents  Comments 


LOCI P+0 
+1 
+ 2 
+ 3 
+4 
+ 5 


I.OCNE+0 
+ ] 
+ 2 

LOCEK+O 


LOC(LOCELM) 

LOCNE 

EH 

EL 

LOSAM 

E1 

e2 

• 

ene 

NE 

BLANK 

HMAXEK 

LCENEK 

UCEKj 

LOCEK 

_ 

1 

: 

LOCEK 

NMAXEK 

LOCNE^IOC  IP  + 5 +NE 


LOCEK“bOCNE+ 3 


NMAXEK  words 


TYPE  6 - continued 


Contents 


LOCNEK+O 


+ 1 


LOCEK^+0 


LOCEK^  +1 


LOCEK.,+0 


+ 1 


LOCEK^+O 


LOCEK  tO 

NMAXEK 


NENEK 

=NE 

NEKX 

1 

nlknmaxek 

X1  at  E1 

X1  at  E2 

Ki  at  enek1 

X2  at  E1 

X2  at  E2 

i 

- 

X2  at  ENEK2 

: 

; 

* 

i 

1 

‘ NMAXEK- 1 
at  E 

niknmaxek-1 

'nmaxek  Jt  El 

• 

' NMAXEK 

| at  enek 

NMAXEK 

> 


NMAXEK  words 


Other  perturbat ions  of  Type  6 
repeat  the  above  layout 


Comments 

LOCNEK=LOCEK+NMAXEK+ 1 
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TYPE  7 - Energy  Distribution  in  Inelastic  Continuum  Scattering 


Contents 


uxnp+o 


+1 


+2 


+ 3 


+4 


+ 5 


ELEN 


1 OCNKP 


LOC (LOCELM) 


NMXENN 


EH 


EL 


NES 


LOSAM 


"NES 


LOCEN 


LOCEN 


LOCEN 


NMXENN 


NEP, 


NEP 


NEP 


NMXENN 


V NMXE! 


NMXENN  words 


J 


NMXENN  words 


Comments 


LLEN=LOCIP+5+NES 


LOCNEP=LLF.N+  NMXENN 


Comments 


TYPE  7 - continued 

Contents 


Other  perturbations  of 
TYPE  7 repeat  the  above 
layout 


TYPE  8 - Cross  Section  for  Inelastic  Level  Excitation 


Contents  Comments 


TYPE  8 - continued 


Contents 


Comments 


LOCMIL 

+LEVMIN 


LEVMIN 


LEVMIN+1 


LEVMIN+IL 


LEV=LEVMIN+IL 

= 1,2,.  (LEVMAX-LEVMIN) 


LEV MAX 


LEVMIN 


LEVMIN 


LEVMIN  at  E, 


millev=0  if  nellev=0 


LEVMIN  at  E 


LEVMIN 


LEVMIN 


LEV  at  E, 


LEV  at  E_ 


LEV  at  E. 


TYPE  8 


- continued 


Contents 


NEL 

LEVMAX 

° LEVMAX 

at 

Ei 

0 LEVMAX 

at 

E2 

° LEVMAX 

at  E 

NEL 

LEVMAX 

Other  perturbations 
of  TYPE  8 repeat  the 
above  layout 


Comments 


TYPE  p - *\ngular  Distribution  of  Discrete  Inelastic  Scattering 
Contents  Comments 


LOCIP+O 

LOC (LOCELM) 

+ 1 

NENS 

+ 2 

EH 

+ 3 

EL 

+ 4 

LKVMIN 

+ 5 

LEVMAX 

+6 

LOCMIL 

+7 

E1 

E2 

• 

enens 

LOCMIL 

H.EVMIN 

mti‘lkvmin 

MIL 


I.EVMIN+  1 


TYPE  9 - continued 


Contents 


U"  MiL+LEV 


mil 


I.KVHIN 


Comments 


LEV=LEVMIN+IL 

IL=  1,2, . . . (LEVMAX-LEVMIN) 


NMU  words  at 


nmu  words  ut  I 


TYPE  10  - Source  Spectrum 


Contents 

Comments 

LOCIP+O 

NES 

1 

+1 

BLANK 

+2 

EH>EHIGH 

+ 3 

EL<ELOW 

+4 

E1 

| 

E2 

* 

LOCIP+3+NES 

enes 
(As/AE) x 

(As/AE)2 

(As/Ae)nes 

(ignore  this 
last  one) 

Other  perturbations 
of  TYPE  10 
repeat  the  above 

format 

End  of  perturbation  data 

BLANK 

LI. NEXT 

BLANK 

Starting  from  location  (LLNEXT-1) 
storage  available  for  next  per- 

turbation data  (e.g.,  sampling 
CHI-  or  ENN-tables ) 
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C.3  - Storage  Description  for  Sampling  CHI-  and  ENN-tables 


of  One  Energy  Band 


Contents 


Comments 


LLNEXT-1 

LLNEXT 

LOCTAX+O 


=LOCSAM 


(if  any) 


NTAE=total  number  of  ENf.-tal . J 

NTAB=total  number  of  CHI-t^; 

LOCTAX=location  of  current 
sampling  table 

LOCTABj  (1=2,  3,...)  = local 

of  next  sampling  table, 


LOCTAX+1 


LOCNIP=location  whert  ti 
number  of  perturbations 
share  the  present  sampling 
resides 


KH= index  of  highest  energy 

g _ unperturbed  energy  table  wher 

1 the  present  sampling  table 

KL-KH+1  begins 

KL=index  of  lowest  energy  of 

unperturbed  energy  table  where 
the  present  sampling  table  ends 


X15  1 


X15  1* 


Xt=0. 


NES^  x 15  words 
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Contents 


Comments 


LOCNIP 
for  table  1 


LOCTAX 
=LOCTAB 1 
LOCTAX+1 
+2 
+ 3 


•LOCNIP 


Kh 


KL 


15 


NIP  words 


15 


’■NIP 

IP 


IP 


NES7=KL-KH+1 


NES?xl5  words 


NIP  words 


IXK'TAX 

J.OCTAB. 


IP 


LOCTAB, 


NEXT 


NIP=total  number  of  IP's  sharing 
the  present  sampling  table 


NOTE:  If  NTAB=0,  it  means  there  is 
no  sampling  x-table  and  the  above 
storage  description  applies  to 
ENN-table  (replace  x by  ENN) . If 
NTAB=1,  then  the  next  table  will  be 
an  ENN-table  if  NTAE^O. 


etc 


APPENDIX  D 


ORGANIZATION  Ot  THE  IPBIN-ARRAY 


The  following  section  describes  the  organization  and  contents  of  the 
IPBIN-array. 


Contents 


Comments 


LOC's  are  pointers  which  g; 
location  which  contains  t:..-  t< 
number  of  perturbations  (!>•*»  i 
affecting  a given  ener  • 
PETAB-array  (cf.  Secti.  i.  .. 

LP  = total  number  of  energy 
in  PETAB-array. 


NIP  words  for  BIN  #1 


17<) 
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BOUNDED  FLUX-AT-A-POINT  (BFAP)  ESTIMATION 
E.I  Introduction 

In  many  instances,  the  availability  of  a point  estimation  capability  re- 
presents a desirable  extension  of  the  ordinary  track  length  or  volume  esti- 
mation method.  However,  unless  special  biasing  procedures  are  employed,  flux- 
at-a-point  estimation  effects  infinite  variance. 

A very  efficient  scheme,  developed  by  Steinberg  and  Kalos * , involves 
biased  selection  of  source  and  collision  points  such  that,  not  only  are  sui  t- 
quelle  l -•  L i.  ht"  estimates  at  point  detectors  made  with  finite  variance,  but. 

estimates  are  made  with  a finite  upper  bound.  This  combination 
of  biased  sel> -it ion  and  "last  flight"  estimation,  the  so-called  bounded  flux- 
at-a-point  (BFAP)  estimation,  is  superior  to  the  so-called  "once-more-collided" 
f lux-at-a-point  (FAP ) estimation^,  because  fewer  histories  are  required  to 
achieve  a given  statistical  accuracy 

The  sections  which  follow  describe  the  mathematical  basis  for,  and  the 
i:.i  implementation  of  the  algorithms  comprising  the  BFAP  estimation  procedvrt 
E : Theory 

. . . 1 . i umt  Estimation 

Tlii  tollowing  procedure  can  be  used  to  make  f lux-at-a-point  estimates.  In 
the  course  of  each  Monte  Carlo  history,  at  the  time  of  position  (source  or 

! i in)  selection,  a virtual  ray  is  traced  from  the  selected  position  to  the 
r.  The  estimator  term,  f,  is  given  by  the  expression  (where  i is  the 
vector  of  the  virtual  ray) : 

2 

f = Wg  ( . ) exp  (~fr>  (s)ds) /4uR 


I HI 


# 

V 

, am 


rl-js.  C il)  1 . X}  i-’A  la 


cLANK-IiOT  ,i'IL  ID 


where  W is  the  particle  weight  at  the  selected  position,  g(fi)  is  the  differ- 
ential directional  distribution  (discussed  below) , a is  the  total  macro- 
\>pic  cross-section,  and  R is  the  distance  between  the  selected  position  and 
the  detector.  The  integration  in  the  exponential  is  along  the  straight  line 
path  from  the  selected  position  to  the  detector. 

The  term  g (ft)  out  of  a source  would  be  constant  for  isotropic  distribu- 
tions; otherwise  it  is  a prescribed  function  of  direction  At  a collision, 
the  situation  rs  more  complicated.  The  collision  procedure  is  usually  carried 
out  in  several  selection  steps:  (1)  position,  (2)  target  element,  (3)  react iui 

process,  (4)  scattering  angle.  The  most  efficient  place  to  carry  out  the  point 
stimation  is  between  steps  (3)  and  (4),  where  W reflects  the  state  of  the 
particle  after  step  (3).  g(ft)  is  then,  simply,  the  (laboratory)  angular  dis- 
tribution of  the  selected  reaction  for  the  selected  target,  and  is  actually  a 

function  of  ft  42,  where  ft  is  the  particle  direction  entering  collision, 
o o 

The  principal  problem  associated  with  such  an  estimation  procedure  is 
2 

i : presence  of  the  R term  in  the  denominator  of  f.  If  the  selected  position 

■m  the  detector,  the  estimator  is  arbitrarily  large  in  value.  As  a means 

. omi'ig  this  problem,  Steinberg  and  Kalos^  devised  a procedure  (discuss,  i 

wiu:re  the  position  selection  is  biased  so  that  W contains  a compensating 
2 

ortional  to  R . This  then  leads  to  an  upper  bound  for  f. 

: . . Bounded  Estimation 

.nee  a full  discussion  of  the  theory  of  bounded  estimation  for  point  de- 
< • is  given  in  Reference  3,  only  the  essentials  will  be  described.  Con- 
: i i ri  around  the  detector  at  x^  (e.g.,  a sphere  of  radius  R^  centered 
di  t.  -tin  , or  a region  of  some  other  shape  containing  sucli  a sphere). 

, position  selection  is  biased  to  be  uniform  in  the  distance  from  the 

i,  t!i.  compensating  weight  factor  contains  a term  proportional  to  R'  , 

. , , ■ k . distance  between  the  detector  and  the  selected  position. 
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When  the  selected  position  x is  a collision  position,  care  must  be  taken 

that  the  distance  to  the  previous  particle  position  x^  is  also  biased  properly 

-2 

to  avoid  a singularity  of  the  form  R^  , where  R^  is  the  distance  between  the 

previous  position  and  the  collision  position.  This  can  be  accomplished  by 

first  biasing  the  direction  out  of  the  previous  collision  so  that  selection 

of  the  angle  a between  (x-x^)  and  (x  -xq)  is  proportional  to  a rather  than  the 

"natural"  cosa  (ignoring  local  anisotropy).  Then  the  position  along  tht  path 

2 

is  selected  from  a probability  density  with  a factor  1/R 

In  Reference  3,  an  implementation  using  a reselection  procedure  was  di 
scribed,  where  x is  reselected  within  the  detector  region  by  reselecting  and 
the  distance  along  the  ray.  (If  the  selected  position  is  a source  position, 

: robl  m does  not  exist,  and  the  method  described  in  Reference  3 applies.) 

The  principal  drawback  of  this  procedure  is  the  need  for  the  availability 
of  the  angular  distribution  data  to  get  the  proper  weight  adjustment  for  the 
selection.  If  the  previous  position  was  also  a collision,  this  angular  dis- 
tribution data  is  a function  of  the  particle  energy  before  the  previous 
collision,  which  might  not  be  available  at  the  time  the  reselection  takes  place, 
because  a different  band  of  cross  section  data  is  in  the  machine.  To  avoid 
this  programming  difficulty,  the  angle  portion  of  the  biasing  is  carried  out 
concurrently  with  the  selection  of  the  direction  out  of  the  previous  position. 
The  biasing  of  the  distance  along  the  ray  is  carried  out  later  in  the  course 
of  selecting  the  collision  position. 

K.  2.3 . Mu  1 t_jj> le  Detectors 

In  in  : n 1 , the  procedures  described  above  can  be  easily  implemented  when 
there  a:-  -v.ral  point  detectors.  The  estimation  procedure  becomes,  essen- 
tially i repetition  of  the  single  detector  scheme,  with  the  added  comi li  ati  n 
of  pos  it,,'  interactions  between  the  detectors.  The  algorithms  utilised  t i 
solve  tin  multi-detector  complication  are  discussed  in  Appendix  F. 


1(3.1 


f 

E ■ 3 Code  Implementation 

The  implementation  of  BFAP  estimation  comprises  several  distinct  aspects: 

(1)  problem  initialization;  (2)  angle  reselection;  (3)  position  selection;  and 
(4)  resolution  of  multiple  detector  conflicts.  The  first  three  aspects  are 
discussed  in  the  sections  that  follow.  Multiple  detector  effects  are  the  sub- 
ject of  Appendix  F. 

E.3.1.  Problem  Initialization 

After  reading  in  user  supplied  information  about  the  point  detectors,  two 
basic  characteristics  are  precomputed  for  each  detector.  The  first  character 
is tic  is  the  so-called  "critical  radius",  which  defines  a sphere  of  influence 

■out  the  detector  point,  also  referred  to  as  the  "critical  sphere".  The  algor- 

★ 

ithm  for  computing  this  radius  is 

; \ = (\)_1 

where  ,f  is  the  total  macroscopic  cross  section  of  the  material  region  designated 

tor  the  k-th  detector,  evaluated  for  a nominal  microscopic  cross  section  of 

1 barn/atom. 

The  second  characteristic  determined  for  each  detector  is  a "critical 

overlap"  flag,  which  is  utilized  in  the  resolution  of  multiple  detector 

' . ff  fx  } are  the  detector  locations,  then  the  alaorithm  for  deter- 

k 

.'t  the  ;c  overlap  flags  (initialized  to  0)  is 

|xk-Xi|2<(Rk+R.)2,  i k;  set  Lk=l.  (E-2) 


■ r may,  optionally,  supply  a value  for  the  critical  radius  directly, 
ribed  in  the  section  on  point  detector  input. 
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E.3.2  Angle  Reselection 

Angle  reselection  is  indicated  whenever  an  originally  selected  (OS) 

ray  intersects  the  critical  sphere  of  a live  detector.  In  addition,  for  post 

scattering  (as  opposed  to  source)  directions,  the  corresponding  post  scat'. r- 

ing  energy  must  be  recomputed.  Furthermore,  for  near  threshold  inelastic  or 

hydrogen  scattering  of  neutrons,  for  which  the  effective  mass  of  recoiling 

nucleus  is  A’<1  [neutron  mass],  a special  treatment  is  required,  since  a rami 

of  scattered  directions  becomes  physically  impossible. 

K.3.2.1.  Angle  Reselection  (General) 

Let  W be  the  originally  selected  (OS)  direction.  The  algoritlim  for  it 

turmi.iing  the  intersection  of  live  critical  spheres  proceeds  as  follows: 

r i h live  detector,  let  W = (x  -x)/|x  -x j , where  x,  is  the  position  of 

a K K K 

the  k-th  detector,  and  x is  the  particle  position.  Let  6^  = W -w  k . Then,  if 
>1  - R^/lx-xJ2,  <Sk>0)  (E— 3 ) 

the  k-th  critical  sphere  is  intersected  by  the  OS  ray.  If  no  live  detector 
satisfies  the  conditions  of  ( E—  3 ) , no  reselection  is  necessary,  and  the  re- 
maining procedure  is  bypassed.  If  more  than  one  live  detector  sphere  is  ir.ter- 
;ected  by  the  OS  ray,  a "multiple  detector  conflict"  (MDC)  exists.  Even  il 
only  one  detector  sphere  is  intersected,  an  MDC  may  have  to  be  resoivi-d  a 
t ■ ' sphere  or  cone  overlap  (see  Appendix  F). 

Ultimately,  subsequent  to  possible  resolution  of  an  original  MDC,  it  m ' 

one  live  detector  sphere,  k,  will  be  intersected  by  the  OS  ray.  The  angle 

* 

i.  selection,  then,  proceeds  as  follows:  choose  an  angle  6 uniformly  in  thi 

■nge  [e  ,0t  1 given  by 

r = min  { R^/ | x-x^ | ,1.0}  (E-4a) 

0fl  = sin  1 (r)  (E-4b) 


For  near  threshold  inelastic  or  hydrogen  scattering,  the  algorithm  of 
(E-4)  is  replaced  by  the  special  procedure  described  in  the  next  section. 

The  reselected  angle  0*  yields  a reselected  (RS)  direction  (or  ra\)  W* 
given  by 


W*  = aw  + gw, 
k 

(E-5a.) 

a = sin  9*/sin  9 

(E-5b) 

g = sin  (e-0*)/sin  0 

-1  - 

(E-5c) 

0 = cos  (W-W^) 

(E-5d) 

To  avoid  the  singularity  when  0=0,  the  following  algorithm  precedes  (£'-5' : 
Let  W = Vi^  = (v^,  v^ , Find  i for  which  jv^j  is  smallest.  Let  new  W be 

the  vector  formed  by  adding  0.01  to  the  i-th  component  of  W , and  normalizing 
to  unity. 

Finally,  the  particle  weight  must  be  adjusted  by  a t i tor 

W = C q(A*)/q(A)  (E-6a) 

a 

wher.  q is  the  sampling  density  for  lab  angle  cosines  ('  and  * are  the  OS 

i:id  RS  values,  respectively),  and 

C = | sin  6*|  (0}l  - 6l)/K0h.9l)  (E-Ob) 

1(9  ,0  ) = COS0  - cos9  , 9 >0 
ML  L H L — 

= 2 - COS0-,  - cos0  ,6  <0  (E-bc) 

H L L 

.(.2.2  Special  Procedure  For  A'^_l  Case 

For  near  threshold  inelastic  and  hydrogen  scattering,  where  thi  (eft  , • r < 

w«  : i!  of  the  recoiling  nucleus  A'^l,  tlie  general  algorithms  (E-3)  and  (F 

i:  modified  as  follows:  in  addition  to  the  constraint  of  (E-3),  the  condition 

W.ft,  < W -ft,  (E-7  ) 

k ok 

i bypass  of  anqular  reselection,  where  W , W,  and  W are  the  pro- 

o k 

’ eui  tli'  direction  to  the  k-th  detector,  respectively. 


If  (E-7)  is  not  satisfied,  the  special  procedure  for  A'<_1  proceeds 


with  the  computation  of  three  angles,  6,6 ,y  in  the  range  (0,it) 


6 = sin  (A1 ) ; 


0<_Q<y/2 


6 = cos  1 (W  *W  ) 

K O 


Y = min  (6,6  ) 


(E-7a) 
(E-7b) 
(E-7c ) 


where  6.,  is  given  by  (E-4b) 

n 


If  6>_y+6  , the  remainder  of  the  special  procedure  is  bypassed.  • •;  - 


wise,  define  b = — , (8>6); 

M 2 


= cos  x j (cos6-cos6cosY)/ (sin6sinY)l  (£<6) 

f “I  J 

Choose  ip  uniformly  in  | 0 , . Compute  a temporary  W,  as  input  tc 


(E-8 ) 


•lecti 


■ ng 


W - A W,  + B W + C (W,  X W ) 
k o k o 


(E-9) 


where 


A = cos6  (1  - cosiM 
B = cos^ 

C = + sin<£ 


(E-lOa) 

(E-lOb) 

(E-lOc) 


The  sign  of  C is  chosen  at  random  with  equal  probability. 

As  a final  preparation  to  the  reselection  procedure,  it  is  necessar 
(•■•fine  " , ti  , and  6.  Substituting  (E-9)  into  (E-5d)  yields  6=5. 


■•  t U = cos^sini 


/ 2„  2 2 

V = j sin  8-sin  t|>sin  6 


2 2 2. 

X = cos  ysin"6+cos  6 


(E-lla) 

(E-llb) 

(E-llc) 


lfi7 


(E-12a) 


Define  0+  in  the  range  (0,2^)  using: 


cos0+  = (cos8cos5+  UV)/X;  ip<  /2 


sin3+  = cosB  - cos<$cos9  + : V<  /2 


(E-i2b) 


sin5cosV 


cos0+  = cos8/cos6  ; ip  = /2 


(E-12c) 


sin0+  = ± J1  ~ 


cos  8+  ; ip 


/2 


(E-12d) 


using  6+  as  defined  above: 


61  = 

max 

(0+,O_) 

; 

(E-l 3a) 

62  = 

min 

(0+,e_) 

; B<6 

(E-l 3b) 

o 

t— ' 

II 

min 

(9+'9_) 

; 8><S 

(E-l 3c) 

II 

CM 

CD 

max 

(e+,0_) 

- 2n  ; B>6 

(E-l 3d) 

9h  " 

min 

( e 1 , y ) 

(E-l 3e) 

9l  = 

max 

(62,-Y) 

(E-l 3f ) 

An  additional 

weight 

adjustment  factor 

C must  be  computed  by : 

it 

M 

iiyv 

c — * 

n-f2- 

ncosB-i^ 

M 

cosy 

“VV 

is  given  by 

(E-6c)  . ft  and  n 

are  angles  between  0 and  " 

n cos  ( (cosy-cos8cos6)/sin6sin5) 
- cos  ^ ( (cos6-cosScosy)/sin8sin'r ) 


(E-14a) 

(E-14b) 
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E.3.3.  Position  Selection 


By  means  of  the  algorithms  (discussed  in  Appendix  F)  designed  to  resolve 
an  MDC,  it  is  possible  for  a particle  ray  to  intersect  at  most  one  sphere  of 
radius  R,  centered  at  a live  detector.  If  a critical  sphere  is  not  inter- 
sected, the  probability  density  for  position  selection  is  given  by  fQ(S) 
for  "ordinary"  tracking,  viz. 


f(S)  = f (S)  = max  u . F . exp  (—  ) 

i 


(E-1S) 


where  i is  the  problem  index  and 

in  = attenuation  coefficient  (for  problem  i) ; 

= problem  dependent  weight; 

S = geometric  distance  along  the  ray. 

In  i ; . t.ice,  the  track  interval  is  subdivided  into  pieces,  with  break  points 
defined  at  points  where  the  index  corresponding  to  the  maximum  (as  computed  in 
E-15  above)  changes. 

If  a live  critical  sphere  of  radius  R is  intersected,  the  probability 

density  f(S)  is  determined  as  follows.  Let  f (S)  be  defined  as  in  (E-15) , 

o 

and  consider  a subinterval  over  which  one  particular  index  i defines  the  maxi- 
mum used  to  obtain  f (S) . Then  f (S)  assumes  one  of  two  forms,  viz. 
o 


f(S)  = Cf  (S) 
o 


(E-16a) 


f(S)  = R f (S  ) /h (S) 
O A 


(E-lLb) 


win  t’  c,  following  algorithms  determine  the  specific  form  for  f(S). 


Let  S.  and  S„  (>S  „)  be  the  end  points  of  the  subinterval  and  let 
ABA  r 


h(S)  = R ~ + (S-S  ) “ 

o o 


(E-16c) 


R distance  between  ray  and  detector; 
o 

distance  along  ray  to  the  point  at  which 


!■  i ■;  the  distance  to  the  deti  itor . 

o 


Noting  that  the  procedure  is  operative  only  when  R <R,  three  breakpoints 
are  defined  initially  (see  Figure  E-l) : 


(E-17a) 

(E-17b) 


(E-l 7c) 


Figure  E.l  - Breakpoint  Definitions 


If  a breakpoint  falls  between  and  S , the  interval  is  divided  at 
the  breakpoint. 

Case  1;  S s S 

L 

In  this  case  the  appropriate  form  for  f(S)  is  given  by  (E-16a) , with  C given 
by  R/Rq,  viz. 


f(S)  = Rf  (s)/R  . 

o o 

Case  2:  S <S<S 

L o 

Compute  r=h(S  )/h(S  ). 

B A 

For  r>l/2 


(E-18) 


f(S>  = R“fo(S)/h(S  ). 


(E-19) 


Otherwise,  for  r<l/2,  further  tests  are  necessary.  If  ^ ..in  . , th  - 
(E-16b)  is  used.  If  the  latter  test  fails,  a new  breakpoint  S is  computed. 


SB  = max  (Sa,Sg) 


(E-20 a) 


where  and  Sg  satisfy 

h<sa)/h(SA)  = 1/2 

WV  = ln  2* 


(E-20b) 


(E-20c) 


Then  f(S)  is  given  by  (E-19)  if  s„  = S ; and  (E-16b)  applies  for  S = S . 

pa  b 

Case  3:  S < S<S 

o H 


The  algorithms  of  Case  2 apply,  with 


r « h(SR)/h(S0) 
C = R2/h(S„) 


(E- 21a) 


(E— 21b) 


and  defined  by 


h(S  )/h(S()  = 0.5. 


(E-21c) 


Wlien  B-Sh,  f(S)  reverts  to  the  ordinary  tracking  density,  f (s)  . In  . 1 1 1 


cases,  the  particle  weight  adjustment  is  given  by 


wi  = UA  exp(-U ,S)/f (s) 


( E - 2 2 ) 


where  i is  t)ie  problem  index. 


APPENDIX  F 


MULTIPLE  DETECTOR  CONFLICTS  (MDC) 


F.l  Introduction 

Since  the  special  biasing  schemes  of  the  BFAP  estimation  procedure  art 
based  on  the  concept  of  detector  spheres  of  influence,  the  simultaneous  pre- 
sence of  several  detectors  defines  potentially  overlapping  regions  of  influ- 
ence. Such  regions  represent  multiple  detector  conflicts  (MDC).  The  MDC  ar 
resolved  via  the  "live  detector"  concept  described  in  the  sections  that  foil 
F.2  The  MDC  Condition 

Relevant  to  the  aspect  of  angle  reselection,  the  MDC  condition  exists 

.r  anv  pa^j . of  live  detectors  satisfies  the  "cone  overlap  test"  (COT) . 

Lor  X,  Xy , and  X , be  the  positions  of  the  particle,  and  the  k-th,  and  the 

1-th  detector,  respectively.  Let  f5.  = W.  'W.,  where  W = (X  — X) / 1 X —X  j , and 

Kic  K 36  K K K 

IV  is  similarly  defined.  Then,  if 


(1 


|x-xkl2 


) (1  - 


|x-V2 


T) 


\R1 


I x-xk ! |x-xs  | 


the  COT  is  satisfied,  and  an  MDC  exists. 

Although  references  to  "sphere  overlaps"  (see  algorithm  (E-2))  . r l 
it  ;>le  intefseetions"  (see  algorithm  ^ E—  3 ) ) also  imnly  the  existence  oi  . 
ui  lb  ’ MDC,  these  tests  comprise  a sub-set  of  the  COT.  Finally,  only  a 
sphere  overlap  is  a significant  MDC  for  position  reselection. 

F 3 Resolution  of  an  MDC 

An  MDC  may  b>  resolved  by  Russian  roulette,  in  which  one  live  detector 
; chosen  for  sampling  purposes.  The  implication  of  this  deactivation  do- 
: on  the  particular  stage  of  the  particle  history  at  which  it  occurs 
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f FA-lr.  tLANK»IiCT  r I L,  L 


F.3.1  Position  Checking  (Source) 


After  the  selection  of  a source  position  at  the  start  of  a history,  the 
resolution  of  an  MDC  is  recorded  by  the  appropriate  setting  of  the  live  de- 
tector "position  flag",  KDLIV.  The  relevent  algorithms  are  best  summarized 
by  a logical  flow  chart: 


f • 

/ INITIAL 
! POSITION 


^ III  SIDE  N 
„ . FECTOR . 


_>■  KDLIV=0 


YES  (IDET) 


YES 

\y 


;dliv=idkt 


continue 


1 

RUSSIAN  ROULETTE 

J PICK  I LET 

■>  KDLIV=-IDET 

j AT  RANDOM 

FIGURE  F-l 

Source  Position  Checking  Algorithms 


Checkin' 


Subsequent  to  position  checking  (both  at  a source  and  a collisi 


point),  there  are  a series  of  direction  checks  aimed  at  resolving  an  a>;<; 


NDLIV 


order  to  distinguish  between  a latent  and  the  current  ray  being  trackt  - 


another  "direction  flag 


IDLIV  assumes  the  role  and  the  value  of  (the  v 


puted)  NDLIV,  whenever  a source  particle  or  latent  is  picked  up  for  • 


As  in  the  setting  of  KDLIV,  the  relevant  algorithms  for  comi 


are  best  summarized  by  a logical  flow  chart 


DLI V 


MATHEMATICAL  APPLICATIONS  GROUP  INC  ELMSFORD  N Y F/G  18/11 

samcep:  a correlated  monte  carlo  neutron  and  gamma  radiation  TR— ETC(U) 

FEB  77  H LICHTENSTEIN.  H STEINBERG.  J BROOKS  DAAD05-75-C-0735 

BRL-CR-330  NL 


3cf  3 

404036469 

a 

END 

DATE 

FILMED 

' 3-77 


AD-A036  469 
UNCLASSIFIED 


By  means  of  the  IDLIV  (NDLIV)  definition,  sampling  along  the  particle 
ray  can  be  affected  by  at  most  one  sphere,  of  radius  centered  at  a detector. 
Thus,  for  IDLIV>0 , detector  sphere  #IDLIV  is  used;  for  IDLIV<0,  # - IDLIV  is 
used  for  sampling  and  a special  weight  adjustment  is  made  (see  Section  F.3.4). 

The  value  of  IDLIV  also  influences  the  position  checking  at  a collision 
point,  as  described  in  the  next  section. 

..3  Position  Checking  (Collision) 

Subsequent  to  the  selection  of  a collision  point  along  the  current  ray, 
which  is  characterized  by  IDLIV  (as  opposed  to  NDLIV,  which  is  computed  .w 
latent),  a value  of  KDLIV  is  computed  for  the  latent  to  be  stored.  This  value 
of  KDLIV  has  the  same  significance  for  the  latent  as  the  KDLIV  which  is  com- 
puted for  the  source  position.  Again,  the  position  checking  algorithms  at  the 
collision  are  best  summarized  by  the  corresponding  flow  charts 
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F.3.4  Interpretation  of  KDLIV 


i 


■ 

i 


As  a result  of  position  checking,  both  for  a source  and  a collision,  the 
setting  of  KDLIV  has  the  following  implications: 

(a)  The  case  KDLIV>0  implies  that  there  was  no  MDC. 

(b)  For  KDLIV<0,  detector  number  (-KDLIV)  is  used  for  sampling 
purposes,  and  a special  weight  adjustment  allows  bounded 
estimation  for  all  detectors'*.  This  weight  adjustment  is 
given  by 

n 

n b,  /Z  b 
k , m 
m=l 

where 

k = | KDLIV | 

bm=  max  (l,Rm2/|x-xJ2) 

(c)  The  value  of  KDLIV  influences  the  direction  checking  preceding  angle 
reselection,  as  described  in  Section  F.3.2  above. 


■I 
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APPENDIX  G 


"LOWER  BOUND"  FOR  BOUNDED  ESTIMATOR 

G.l  Introduction 

Steinberg^  introduced  a collision  selection-reselection  technique  which 
guarantees  that  the  estimates  to  a point  detector  are  bounded  from  above.  If 
optimum  importance  sampling  is  used,  all  estimates  to  the  detector  should  be 
of  the  order  of 


Q = W(D,E)/4ttR 


(G-l) 


where  W(D,E)  is  the  region  and  energy  dependent  weight  in  the  detector  i go  i , 

at  the  energy  of  the  estimate,  and  Rc  is  the  radius  of  the  "critical  sphere" 

★ 

around  the  detector  (distance  corresponding  to  one  m f.p.)  . 

As  the  importance  sampling  actually  used  is  often  quite  different  from 
the  optimum  one,  the  individual  estimates  will  fluctuate  around  the  mean.  In 
order  to  speed  up  the  calculation,  we  introduce  a low  value  cutoff,  below 
which  estimates  are  Russian-rouletted.  The  low  value  cutoff  has  been  chosen  to 
be  QxFq,  where  Fq  is  the  same  as  F^,  an  input  cutoff  value  relative  to  unity, 
used  in  other  parts  of  the  code. 

G 2 Implementation 


The  quantity  to  be  scored  is  of  the  form: 
S = Ge 


(G-2) 


where  G reflects  all  biasing  and  the  1/r  factor,  and  is  the  optical  dis- 
tance from  a point  to  the  detector,.  We  rewrite  (G-2)  as 
-(X+X  ) 

S = Qe  o (G-3) 

where  Q is  given  by  Eq.  (G-l)  and 

Xq  = -log  (G/Q)  (G-4) 


To  safeguard  aqainst  poor  importance  sampling,  R = max  (R  ,r), 

c c 

where  r is  the  distance  between  particle  position  and  detector. 
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We  further  define  a m.f.p.  cutoff  value 


X = -log  (F  ) 
c o 

(G-5) 

and  replace  the  test 

S < QxFq 

(G-6) 

by  the  test 

>* 

+ 

0 

V 

>* 

o 

(G-7) 

The  optical  distance  X is  calculated  by  tracking  as  a cumulative  sum  of  non- 
negative terms. 

X = Xn,  where 

k k 

XK  = I X. 

i 


with  X.  > 0 for  i > 0. 

i 

Tlie  test  (G-7)  is  performed  by  testing 


for  k=0,l... 


(G-8) 


until  either  the  test  is  satisfied,  or  k=n,  whichever  occurs  first.  If  the 

test  is  satisfied,  a game  of  chance  is  played.  With  probability  1-F  , the 

o 

estimator  becomes  0 and  the  tracking  is  terminated.  With  probability  Fq,  the 

k k 

estimator  is  multiplied  by  1/F  (which  is  achieved  by  replacing  X by  X -X  ), 

o c 

the  tracking  and  test  (G-8)  continues  for  succeeding  values  of  k. 


200 


APPENDIX  H 


ENERGY  RANGE  HIERARCHY 

The  principal  energy  range  for  a given  SAMCEP  run  is  the  Monte  Carlo 
tracking  range,  as  specified  in  the  SAMCAR  input  by  EHIGH  and  ECUT  (see 
Item  11,  Section  5.3).  The  other  major  energy  ranges  are  defined,  relative 
to  (EHIGH,  ECUT),  as  follows: 

MONTE  CARLO 

a)  EOUT,  > EHIGH  > ECUT  > EOUT 

1 LAST 

where  (EOUT^EOUT^^)  is  the  output  energy 
range  for  scoring  (Item  12,  Section  5.3); 

b)  ES  > EHIGH  > ECUT  > ES 

-L  LAST 

where  (ESjyES^  ) is  the  source  energy  range 
(Item  25,  Section  5.3); 

c)  EWTAB,  > EHIGH  > ECUT  > EWTAB. 

1 LAST 

where  (EWTAB^EWTAB^^)  is  the  energy  importance 
range  (Item  17,  Section  5.3); 

NOTE:  There  is  no  hierarchy  established  among 

EOUT,  ES,  and  EWTAB. 

DATA 

d)  EEHIGH  > EHIGH  > ECUT  > EELOW 

where  (EEHIGH, EELOW)  is  the  maximum  energy  range 
defined  by  perturbation  input  (Item  1,  Section  3.2); 

e)  EX  > EEHIGH  > EELOW  > EX 

1 LAST 

where  (EX^  'EXLAST)  is  the  minimum  cross  section 
energy  range  (SAM-X  output)  utilized  in  the  problem 
(including  type  2 perturbations). 
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APPENDIX  I 


OPERATING  INSTRUCTIONS  FOR  SAMCEP 


SAMCEP  is  currently  operational  under  the  SCOPE  2.1  system  on  the  CDC 
7600.  A card  deck  configuration,  assuming  the  five  executable  programs 
comprising  the  SAMCEP  system  reside  on  a permanent  file  SAMCEP,  is  shown  belov.: 
CARD  CONTENTS  COMMENTS 

SBLMJ, STMFZ.  7600  JOB  card. 


ATTACH (SAMIN, SAMCEP, CY=1) 
SAMIN. 

REWIND (TAPE8) 

ATTACH ( TAPE 1 1 , NEDT , C Y= 1 ) 
ATTACH ( SAMSAM , SAMCEP , CY=  2 ) 
SAMS AM. 

REWI ND ( TAPE  8 , TAPE 1 0 , TAPE 1 2 ) 


RETURN (TAPE 11 ) 

ATTACH ( S AMCAR, SAMCEP, CY= 3) 
SAMCAR  . 

RETURN (TAPE10) 

RETURN (TAPE12) 


Execute  the  neutron  perturbation  proc^sst 
Processed  perturbation  data  tape. 

Neutron  element  data  generated  by  SAM-X. 

Execute  the  neutron  transport  preprocessor. 
TAPE8  is  the  modified  perturbation  data  ta; 
TAPE10  is  the  organized  data  tape  (ODT  ou_pa 
of  BAND);  TAPE12  is  the  sampling  CHI-and/or 
ENN-table  tape  (if  any) . 


Execute  the  correlated  Monte  Carlo  transit*- 
code  for  the  primary  neutron  problem. 

The  neutron  ODT  may  also  be  saved  for  futur 
neutron  runs. 

The  sampling  tape  may  also  be  saved  for 
future  neutron  runs. 
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L 


J ^1<jT 


CARD  CONTENTS 


COMMENTS 


REWIND (TAPE8 , TAPE14 , TAPE16) 


A I i’;V  I ( i AMOOT , SAMCEP  , CY=4 ) 
■>AMC  JT  (TAPEl 6) 

RETURN (TAPE16) 

RETURN (TAPE4) 


ATTACH ( TAPE1 1 , GEDT , CY= 1 ) 
ATTACH (TAPEl 2 , GPDT, CY=1 ) 
ATTACH ( SAMGAM , SAMCEP , CY=  5 ) 

SAMGAM. 

REWIND (TAPES, fAPE10,TAPEl5) 


RETURN (TAPEl 2 ,TAPE14) 


TAPE8  is  needed  for  a secondary  gamma 
problem;  TAPE14  may  contain  non-elastic 
neutron  interactions  from  which  a 
secondary  gamma  source  is  generated;  TAPE 16 
is  the  statistical  aggregate  tape  for  the 
primary  neutron  problem. 

Execute  the  edit  code  (TAPE 16  is  denoted 
TAPEl  internally). 

TAPE4,  the  statistical  tape  generated  by 
SAMOUT,  may  be  saved  for  future  reruns  of 
SAMOJT. 

Gamma  element  data  tape  generated  by  SAM-X. 
Gamma  production  data  tape  generated  by  SAM-X. 

Execute  the  secondary  gamma  preprocessor. 

TAPER  has  been  modified  by  SAMGAM,  retaining 
only  concentration  perturbations  from  the 
primary  neutron  problem;  TAPE10  is  now  the 
organized  gamma  cross  section  data;  TAPE15 
is  the  external  secondary  gamma  source  tape 
generated  by  SAMGAM  from  the  TAPEl 4,  TAPEl 2, 
and  user  supplied  data. 


i 
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CARD  CONTENTS 


COMMENTS 


SAMCAR. 

REWIND (TAPE16) 

SAMOUT (TAPE16) 

(eor) 

data  deck  for  SAM IN 
(eor) 

data  deck  for  SAMSAM 
(eor) 

data  deck  for  primary  neutron  SAMCAR 
(eor) 

data  deck  for  primary  neutron  SAMOUT 

(eor) 

data  deck  for  SAMGAM 
(eor) 

data  deck  for  secondary  gamma  SAMCAR 
(eor) 

data  deck  for  secondary  gamma  SAMOUT 
(eof ) 

- End  of  Job  Deck  - 


Execute  the  correlated  Monte  Carlo  transport 
code  for  the  secondary  gamma  problem. 

TAPE 16  now  contains  statistical  aggregates 
for  the  secondary  gamma  problem. 

Edit  the  secondary  gamma  results. 
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These  programs  can  be  run  in  separate  jobs  by  saving  the  appropriate  tape 


files.  The  following  table  specifies 
execution: 

TO  START  WITH 

SAMIN 

SAMSAM 

SAMCAR 


SAMOUT 

SAMGAM 


the  required  files  for  each  program 

FILES  REQUIRED 
None 

TAPE 8 from  SAMIN;  TAPE11  (NEDT) 
from  SAM-X. 

TAPE8,  TAPE10,  TAPE12  from  SAMSAM 
(for  primary  neutron); 

TAPE8,  TAPE10,  TAPE15  from  SAMGAM 
(for  secondary  gamma)  . 

TAPE 16  from  SAMCAR  or  TAPE 4 from 
previous  SAMOUT. 

TAPE 8 from  SAMSAM;  TAPE14  from  primary 
neutron  SAMCAR;  TAPE11  (GEDT) , TAPE 12 
(GPDT)  from  SAM-X. 
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Philadelphia,  PA  19137 

2 Commander 

US  Army  Picatinny  Arsenal 
ATTN:  SARPA-ND-C-T 

JAWT1P/SARPA-RT-S 
Dover,  NJ  07801 

1 Commander 

US  Army  Nuclear  Agency 
Fort  Bliss,  TX  79916 

1 Commander 

US  Army  Harry  Diamond  Labs 
ATTN:  DRXDO-TI 

2800  Powder  Mill  Road 
Adelphi , MD  20783 

1 Director 

US  Army  TRADOC  Systems 
Analysis  Activity 
ATTN:  ATAA-SA 

White  Sands  Missile  Range 
NM  88002 

1 AFWL  (SUL) 

Kirt land  AFB , NM  87117 
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1 Director 

Lawrence  Livermore  Laboratory 
ATTN:  Tech  Info  Div 

P.  0.  Box  808 
Livermore,  CA  94S50 

1 Director 

Los  Alamos  Scientific  Lab 
ATTN:  Repts  Lib 

P.  0.  Box  1663 
Los  Alamos,  NM  87544 

3 Mathematical  Applications 
Group,  Inc. 

3 Westchester  Plaza 
Elmsford,  NY  10523 

1 Science  Applications,  Inc. 
ATTN:  Dr.  E.  Straker 

1200  Prospect  Street 
La  Jolla,  CA  92037 


Aberdeen  Proving  Ground 


Marine  Corps  Ln  Ofc 
Dir,  USAMSAA 
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